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SUMMARY  OF  RESULTS 

A  dipole  inversion  method  was  developed  and  successfully  applied  to  the  discrimination 
and  identification  of  unexploded  ordnance  using  total-field  magnetometry.  A  black- white 
classification  of  anomalies  as  ordnance/non-ordnance  is  not  possible  using  magnetics 
because  some  non-ordnance  items  are  indistinguishable  from  UXO.  We  found  that  the 
best  strategy  is  to  rank  items  according  to  the  likelihood  they  are  UXO.  This  ranking 
can  be  achieved  by  several  methods;  the  best  was  to  rank  on  the  basis  of  the  amount 
of  remnant  magnetization  required  to  make  the  anomaly  match  one  of  the  ordnance 
items  in  a  predefined  library.  The  discrimination  method  had  the  potential  to  reduce 
the  number  of  excavations  at  an  impact  site  in  Montana  (Guthrie  Road)  by  half,  yet 
still  yielded  all  the  ordnance.  The  method  requires  complete  or  at  least  partial  shock 
demagnetization  of  ordnance,  otherwise  it  is  not  possible  to  reliably  rank  the  anomalies. 

Error  analysis  indicates  that  recovered  dipole  moments  are  relatively  well  constrained. 
However,  identification  of  the  anomaly  source  is  difficult  using  the  recovered  dipole  mo¬ 
ment  due  to  non-uniqueness.  The  dipole  moment  is  the  product  of  magnetization  with 
volume  and  a  change  in  volume  can  be  compensated  by  a  change  in  magnetization 
by  varying  the  orientation  of  the  item  relative  to  the  Earth's  field.  By  recognizing  that 
there  are  a  finite  number  of  ordnance  types,  correct  identification  at  Guthrie  Road  could 
be  achieved  about  55%  of  the  time.  Attempts  to  constrain  the  orientation  (and  hence 
ordnance  dimension)  by  recovering  higher  order  moments  (quadrupole  and  octupole)  of 
the  body  did  not  improve  our  ability  to  discriminate  or  identify. 
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1  INTRODUCTION 


Unexploded  ordnance  pose  a  significant  public  safety  hazard  in  many  parts  of  the  world. 
They  occur  on  or  near  the  surface  and  down  to  depths  of  several  meters  below  the 
ground.  One  source  of  UXO’s  are  armed  conflicts  both  old  and  more  recent  in  regions 
such  as  Europe  (World  Wars  I  and  II),  South-East  Asia,  Afghanistan,  Yugoslavia  and 
parts  of  Africa.  They  are  also  a  significant  problem  in  countries  such  as  Canada  and 
the  United  States,  where  they  are  present  in  areas  used  for  military  training  and  firing 
ranges;  largely  as  vestiges  of  World  War  II  and  the  Cold  War.  It  is  estimated  that  15 
million  acres  of  the  United  States  are  contaminated  by  UXO  (FAC,  1996),  with  a  clean¬ 
up  time-frame  in  the  decades  and  a  staggering  cost  estimate,  using  existing  technology, 
in  the  tens  to  hundreds  of  billions  of  dollars  (GAO,  2001). 

The  most  well  established  techniques  for  ordnance  detection  are  magnetics  and  elec¬ 
tromagnetics  (Bulter  et  al.,  1998).  These  methods  are  very  effective  at  locating  buried 
metallic  objects  such  as  UXO.  However,  discriminating  between  intact  UXOs  and  unin¬ 
teresting  objects  such  as  metallic  debris  and  shrapnel  is  a  significant  challenge  (Bulter 
et  al.,  1998).  Often  little  or  no  effort  is  invested  in  discrimination  and,  consequently, 
many  holes  are  excavated  for  each  ordnance  item  recovered.  For  example,  Putnam 
(2001)  reports  that  from  49,521  anomalies  excavated  at  Kaho’olawe  only  3%  (or  1,486) 
turned  out  to  be  UXO.  In  this  report  we  analyze  magnetics  data  collected  at  Guthrie 
Road,  Helana  Valley,  Montana  (Youmans  and  Daehn,  1999).  840  holes  were  excavated 
with  83  UXOs  recovered,  thus  there  was  a  10%  success  rate.  It  is  perhaps  a  reflection 
of  the  state  of  our  ability  to  discriminate  that  this  survey  was  considered  a  success  by 
all  parties  involved. 

The  simplest  form  of  discrimination  used  is  to  threshold  the  collected  data  in  some 
way  (for  magnetics  this  is  usually  achieved  by  first  calculating  the  analytic  signal  and 
then  thresholding  Roest  et  al.  (2001)).  As  the  threshold  level  is  decreased,  the  number  of 
items  that  need  to  be  excavated  increases.  As  demonstrated  by  the  Ordnance  Detection 
and  Discrimination  Study  at  the  Former  Fort  Ord  (Asch  and  Staes,  2001),  for  such 
techniques  to  yield  a  high  probability  of  detection,  a  very  high  false  alarm  rate  results. 

Location  of  UXOs  involves  the  successful  detection  and  identification  of  compact 
metallic  objects.  Therefore,  a  better  method  for  discrimination  is  to  parameterize  the 
characteristics  of  the  target  in  some  way  (e.g.  by  its  location,  orientation,  shape,  size 
and  material  properties)  and  develop  a  good  forward  model  that  allows  predicted  data 
to  be  generated  from  a  given  target  realization.  Inversion  can  then  be  used  to  find  the 
set  of  target  parameters  that  best  fits  the  data.  The  target  parameters  are  then  used 
to  make  decisions  regarding  discrimination  (UXO/non-UXO)  and  identification  (what 
is  the  ordnance  type). 

A  previous  report  (Pasion  and  Oldenburg,  2001)  considered  the  development  of  an 
inversion  methodology  for  time-domain  electromagnetic  data.  In  this  report,  we  con¬ 
sider  the  development  of  a  methodology  for  total-field  magnetics  data.  Previous  work 
in  this  area  has  been  conducted  by  (amongst  others)  the  Canadian  Defence  Research 
Establishment  (McFee,  1989)  and  the  Naval  Research  Laboratory  (Nelson  et  al.,  1998) 
and  private  companies. 

Spheroids  have  been  proposed  as  an  approximate  parameterization  of  ordnance  by 
several  authors  (McFee,  1989;  Altshuler,  1996;  Bulter  et  al,  1998).  While  the  spheroid 
does  not  capture  the  top-bottom  asymmetry  of  many  ordnance  items,  close  agreement 
between  observed  anomalies  over  test  stands  and  spheroid  fits  have  been  demonstrated 
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(McFee,  1989;  Bulter  et  al.,  1998).  Furthermore,  the  magnetic  anomaly  from  a  solid 
spheroid  has  been  shown  to  be  very  similar  to  a  hollow  spheroid  (Altshuler,  1996). 
Therefore,  we  will  use  solid  spheroids  as  the  basis  of  our  magnetic  modelling  of  UXO’s. 

The  response  of  a  spheroid  (in  fact  any  compact  body)  can  be  decomposed  into  a 
series  of  moments  by  a  multipole  expansion  (Stratton,  1941).  The  response  of  the  dipole 
component  dominates  the  anomalies  caused  by  most  buried  ordnance  due  to  the  rapid 
falloff  with  distance  of  the  other  components.  Consequently,  the  main  inversion  routine 
developed  in  this  report  attempts  to  recover  the  dipole  moment  and  position  that  best 
matches  an  observed  anomaly.  The  recovered  dipole  moment  in  then  used  to  rank  items 
according  to  the  likelihood  they  are  UXO.  We  prefer  such  a  ranking  over  a  definitive 
statement  such  as  ordnance  or  non-ordnance  because  it  is  not  possible  to  make  such  a 
black  and  white  distinction. 

The  recovered  dipole  moment  can  also  be  used  for  ordnance  identification.  There  is 
some  ambiguity  in  this  process  because  the  dipole  moment  is  the  product  of  magneti¬ 
zation  with  volume.  Due  to  self-demagnetization  effects  the  magnetization  can  change 
significantly  as  the  orientation  of  the  body  relative  to  the  Earth’s  field  changes.  Thus 
orientation  can  be  traded  with  volume  to  produce  the  same  or  similar  dipole  moment. 
By  recognizing  that  there  are  a  finite  number  of  ordnance  types  present  in  an  area  (we 
build  up  a  library),  this  ambiguity  can  be  reduced  significantly. 

To  attempt  to  constrain  the  orientation  we  develop  two  further  inversion  routines. 
The  first  uses  the  octupole  component  of  the  spheroid,  and  finds  the  optimum  orientation 
and  fit  for  each  of  the  items  in  our  library.  That  item  with  the  best  fit  is  assumed  to  be 
the  source  of  the  anomaly.  The  second  method  aims  to  recover  the  dipole  and  quadrupole 
moments  of  an  arbitrary  axially  symmetric  body. 

No  inversion  methodology  would  be  complete  without  an  analysis  of  the  uncertainties 
in  the  recovered  parameters.  The  last  section  of  the  report  considers  the  development  of 
a  linearized  error  analysis.  We  also  investigate  the  distribution  of  the  residuals  so  that 
we  can  test  our  assumption  of  Gaussian  statistics. 

2  MAGNETIC  FIELD  OF  A  SPHEROID 

The  forward  modelling  task  may  be  stated  as  follows;  calculate  the  magnetic  anomaly  at 
the  arbitrary  point  (x,y,z)  given  a  spheroid  of  diameter  a,  eccentricity  e  (i.e.  length  = 
ae)  and  permeability  fi  oriented  at  an  angle  of  8  degrees  to  the  horizontal  with  an 
azimuth  of  <f>  degrees  clockwise  from  North,  at  a  depth  z  below  the  surface  and  at 
the  horizontal  position  (x,  y).  The  geometry  of  this  situation  is  shown  in  Figure  1. 
Full  characterization  of  the  anomaly  caused  by  a  spheroid  requires  the  Earth’s  field 
plus  the  eight  parameters  (x,  y,  z,  0,  <£,  a,  e,  y).  For  the  Earth’s  field  we  use  geographical 
coordinates;  B0  =  {Box,  Boy,  Boz),  i.e.  x  is  positive  to  the  East,  y  is  positive  to  the 
North  and  z  is  positive  upwards.  Note  that  this  differs  from  the  definition  used  to 
specify  the  International  Geomagnetic  Reference  Field  (IGRF)  where  x  is  positive  to 
the  North,  y  is  positive  to  the  East  and  z  is  positive  downwards. 

The  magnetic  anomaly  of  a  buried  ferrous  item  arises  from  both  remnant  Mrem  and 
induced  magnetization  Min;  the  total  magnetization  M  is  then 

NI  —  M rern  "h  (1) 

Remnant  magnetism  is  present  even  in  the  absence  of  an  inducing  field  and  is  due  to 
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Figure  1:  Geometry  for  modelling  the  response  of  a  spheroid. 


magnetic  moments  being  looked  into  alignment  once  the  steel  casing  cools  to  below  the 
Curie  temperature  in  the  presence  of  an  external  field.  We  will  discuss  this  type  of 
magnetism  latter  but  for  now  we  will  concentrate  on  the  induced  magnetism  that  arises 
because  magnetic  domains  in  a  ferrous  material  tend  to  align  with  the  direction  of  the 
Earth’s  field.  The  ease  with  which  the  moments  align,  and  hence  the  strength  of  the 
magnetism,  depends  on  the  magnetic  permeability  of  the  steel.  For  a  compact  body 
such  as  a  spheroid,  demagnetization  effects  become  important.  This  phenomena  refers 
to  the  extent  that  the  induced  field  is  reduced  due  to  the  shape  of  the  spheroid.  It 
arises  due  to  the  boundary  conditions  that  the  field  must  satisfy  across  a  discontinuity 
in  magnetic  permeability.  A  solution  of  the  boundary  value  problem  (Stratton,  1941) 
shows  that  the  demagnetization  factors  are 


F  =  —  - 

1  1  +  —  l)/2  (Z) 

where  fir  is  the  relative  permeability  (/i  =  /nr/t0,  where  /j,0  =  in  x  10-7  H/in  is  the 
permeability  of  free  space)  and  the  cv,’s  are  demagnetization  factors  which  depend  on 
the  shape  of  the  spheroid  (see  Appendix  A).  Once  the  relative  permeability  exceeds  a 
few  hundred  units,  the  induced  magnetization  becomes  virtually  independent  of  suscep¬ 
tibility.  This  can  easily  be  seen  from  Equation  2  with  large  p,r  because  then 


This  allows  us  to  eliminate  fir  from  our  calculations  because  steel  typically  has  suscep¬ 
tibilities  of  several  hundred  to  over  a  thousand  (Lide,  2001). 

The  demagnetization  factors,  together  with  the  Earth’s  field  B0  (in  spheroid  centered 
coordinates),  determine  the  strength  of  the  induced  magnetization, 


Min  =  Vo  !FB0 


(4) 
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where  we  have  written  F  as  a  3  x  3  diagonal  matrix  with  (F2,  F2,  F3)  along  the  diagonal 
(remember  F\  —  F2).  An  important  consequence  of  demagnetization  is  that  the  induced 
magnetism  can  change  significantly  with  orientation.  For  instance,  with  e  =  4  we  find 
F2  =  2.1  and  F3  =  12.7  so  that  the  magnetism  when  the  spheroid  semi-major  axis  is 
aligned  with  the  field  is  around  6  times  greater  than  when  it  is  perpendicular. 

Equation  4  returns  Min  in  spheroid  centered  coordinates  and  assumes  B0  is  also 
spheroid  centered.  To  use  geographical  coordinates  we  utilize  the  Euler  rotation  tensor 


A  = 


COSXp 

sin  9  sin  xp 
cos  6  sin  xp 


—  sin  xp 
sin  9  cos  xp 
cos  9  cos  xp 


0 

COS0 

—  sin# 


(5) 


to  rotate  B0  to  spheroid  centered  coordinates  and  then  use  its  inverse  Ar  to  rotate  Min 
to  geographical  coordinates, 


Min  =  Mo  'A^ABo  (6) 

With  this  uniform  magnetism  the  magnetic  field  of  the  spheroid  is  given  by  the  expression 

B  =  <7> 

where  r'  is  the  distance  from  the  source  to  the  observation  point  (Figure  1).  An  exact 
solution  of  this  equation  using  prolate  spheroidal  harmonics  is  known  (Stratton,  1941) 
but  we  chose  to  use  a  multipole  expansion;  the  details  can  be  found  in  Appendix  A  and 
McFee  (1989).  The  first  non-zero  moment  is  the  dipole  m  which  is  given  by 

m  =  VMin  =  — e<z  (8) 

where  V  is  the  volume  of  the  spheroid.  The  magnetic  field  from  the  dipole  term  is  then 

B  =  4^(^[x-m]x-m)  (9) 

The  octupole  is  the  next  non-zero  moment  and  is  a  rank  3  tensor  with  27  components. 
Symmetry  reduces  the  number  of  independent  components  to  6.  The  specification  of 
the  octupole  term  is  quite  complicated,  so  we  leave  the  details  to  Appendix  A  and  just 
note  that  its  field  dies  off  as  r“5. 

The  rapid  falloff  of  the  octupole  term  with  distance  means  that  once  the  sensor 
distance  exceeds  a  few  body  lengths,  the  field  is  essentially  dipolar.  For  example,  Fig¬ 
ure  2a  plots  the  dipole  and  octupole  fields  from  a  spheroid  model  of  a  105  mm  ordnance 
(a  =  0.105m,  e  =  3.7)  orientated  at  45°  to  the  vertical,  along  horizontal  lines  at  a  height 
of  ae  and  2ae  above  the  ordnance.  The  responses  are  normalized  by  the  maximum  dipole 
response  at  that  height.  At  a  height  of  ae  the  peak  of  the  octupole  term  is  a  little  over 
15%  of  the  peak  of  the  dipole  term.  For  the  height  2ae  the  octupole  peak  is  less  than 
5%  the  size  of  the  dipole  peak.  Figure  2b  shows  how  this  ratio  rapidly  decreases  with 
increasing  height  above  the  ordnance. 

Considering  that  real  magnetic  profiles  are  contaminated  by  noise,  observing  the 
octupole  response  in  the  presence  of  the  dipole  response  will  generally  be  very  difficult. 
It  may  be  possible  when  the  sensor  is  very  close  to  the  ordnance,  or  when  the  noise 
levels  are  low  and  the  spatial  positioning  extremely  accurate.  We  now  examine  the 
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Figure  2:  (a)  Dipole  and  octupole  terms  at  one  and  two  body  lengths  distance;  (b)  Logio 
of  the  ratio  of  the  octupole  to  dipole  terms  for  a  105  mm  projectile  at  45°  inclination 


consequence  of  the  dipolar  nature  of  the  observed  fields  on  both  discrimination  and 
identification. 


2.1  Discrimination  using  magnetometry 

The  dipole  moment  of  a  sphere  is  constrained  to  lie  in  the  direction  of  the  applied  field. 
On  the  other  hand,  the  angle  the  induced  dipole  moment  of  a  spheroid  makes  with  the 
Earth’s  field  will  vary  with  its  orientation.  We  show  in  Appendix  B  that  the  angle  ^ 


5 


between  the  induced  moment  and  the  applied  field  will  be 


ip  =  arccos 


F2  sin2  8  4-  Fs  cos2  8 
y/  F$  sin2  0  4-  F$  cos2  6 


(10) 


where  8  is  the  angle  of  the  spheroid  semi-major  axis  with  the  Earth’s  field  and  F2  and 
Fs  are  the  demagnetization  factors  of  the  spheroid.  In  Figure  3a  we  show  the  variation 
in  the  angle  ip  as  the  orientation  8  changes  for  a  spheroid  with  an  aspect-ratio  of  5. 
The  maximum  angle  is  bounded  and  dependent  on  the  shape  of  the  spheroid  (Altshuler, 
1996).  The  spheroid  orientation  8  that  makes  the  largest  angle  (see  Appendix  B)  is 


8  —  arctan 


(  lFj  +  FzF$-2F^f\ 
Fi  +  F2Fi-2F3FZ) 


(11) 


The  maximum  angle  for  a  given  aspect  ratio  can  then  obtained  by  substituting  Equa¬ 
tion  (11)  into  (10).  The  variation  of  this  maximum  angle  with  shape  is  shown  in  Fig¬ 
ure  3b.  Most  UXO’s  have  eccentricities  between  3  and  7,  so  that  55°  is  an  upper  bound 
on  the  angle.  This  suggests  a  possible  method  for  discrimination  (Altshuler,  1996). 

When  an  ordnance  hits  the  ground  but  does  not  explode,  the  impact  is  usually  suf¬ 
ficient  to  cause  shock  demagnetization;  the  shock  of  impact  causes  magnetic  domains 
within  the  ordnance  to  become  randomly  aligned,  thus  erasing  any  remnant  magnetiza¬ 
tion.  On  the  other  hand,  when  the  ordnance  explodes,  the  shrapnel  is  heated,  deformed 
and  cooled  in  the  presence  of  the  Earth’s  field  and  can  acquire  permanent  magnetism. 
Putting  all  of  the  above  together  it  is  conjectured  that  ordnance  should  not  have  a  dipole 
moment  that  exceeds  an  angle  of  60°  from  the  Earth’s  field.  The  success  of  this  method 
of  discrimination  in  reducing  false  alarms  has  been  demonstrated  in  Nelson  et  al.  (1998). 
However,  they  did  note  that  subsequent  events  (such  as  lightning  strikes)  can  cause  the 
ordnance  to  acquire  a  remnant  magnetization  and  hence  degrade  the  performance  of 
the  discrimination  method.  Also,  the  assumption  that  ordnance  are  fully  demagnetized 
on  impact  may  be  violated  in  practice.  We  shall  return  to  this  point  later  but  for  now 
assume  that  our  initial  conjecture  is  valid. 


2.2  Ordnance  identification  using  magnetometry 

We  now  turn  to  the  issue  of  identifying  ordnance  items  from  their  dipole  field.  For  a 
given  spheroid,  the  induced  dipole  moment  will  be  dependent  on  the  angle  8  that  the 
spheroid  axis  makes  with  the  Earth’s  field.  There  is  no  azimuthal  dependence  due  to  the 
spheroid’s  symmetry.  By  varying  the  orientation  8  we  are  able  to  trace  the  change  in 
dipole  magnitude  and  angle  with  respect  to  the  applied  field.  Figure  4  shows  these  curves 
as  polar  plots  for  five  different  ordnance  items;  a  60  mm  mortar,  a  76  mm  projectile, 
an  81  mm  mortar,  a  105  mm  projectile  and  a  155  mm  projectile.  Fitting  a  dipole  to 
observed  data  will  produce  a  single  point  in  this  polar  plot  that  needs  to  be  matched  to 
an  ordnance  item.  The  curves  for  the  76  mm  projectile  and  the  81  mm  mortar  are  very 
similar  in  certain  parts  of  the  plot,  indicating  that  discriminating  between  these  items 
may  be  very  difficult.  Further,  certain  orientations  lead  to  identical  dipole  moments  of 
both  these  items  to  the  60  mm  mortar  and  the  105  mm  projectile. 

The  analysis  in  the  last  paragraph  indicates  that  there  can  be  ambiguity  in  identi¬ 
fication  using  magnetometry.  This  fact  is  rigorously  established  in  Appendix  C  where 
it  is  shown  that  for  a  given  spheroid  at  a  particular  orientation,  there  are  an  infinite 
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Figure  3:  (a)  Change  in  the  moment  angle  as  the  orientation  varies  for  an  aspect  ratio 
of  5;  and  (b)  Maximum  angle  between  the  induced  dipole  moment  of  a  spheroid  and  the 
Earth’s  field 

number  of  other  spheroids  that  could  have  produced  the  same  dipole  moment.  For  ex¬ 
ample,  Figure  5  and  Table  1  show  the  family  of  spheroids  that  can  produce  the  same 
dipole  moment  as  a  105  mm  projectile  orientated  at  45°  to  the  Earth’s  field.  The  90 
mm  projectile  lies  very  close  to  this  curve  emphasizing  the  difficulty  in  separating  these 
two  items.  The  folding  over  of  this  curve  to  encompass  spheroids  of  larger  dimension  is 
of  more  concern  to  classification;  the  155  mm  lies  quite  close  to  this  branch  of  the  curve. 
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Figure  4:  Magnitude  of  dipole  moment  (Am2)  and  angle  from  Earth’s  field  for  60,  76, 
81,  105  and  155  mm  projectiles 

There  is  also  a  family  of  oblate  spheroids  that  could  produce  the  same  dipole  moment  j 
as  could  objects  of  different  shapes.  This  inherent  ambiguity  has  a  significant  impact 
on  the  methods  we  develop  for  ordnance  classification  in  the  remainder  of  this  report. 

3  INVERSION  METHOD  FOR  PARAMETER  ESTIMA¬ 
TION 

The  previous  section  demonstrated  that,  in  many  situations,  it  would  be  difficult  to 
recover  more  than  the  dipole  component  of  a  spheroid’s  magnetic  anomaly.  Therefore, 
rather  than  attempting  a  direct  inversion  for  the  spheroid  location  and  dimension,  we 
apply  a  two-step  procedure: 

1.  Invert  for  the  best-fitting  dipole  moment  and  location; 

2.  Use  the  recovered  dipole  moment  for  discrimination  and  ordnance  identification. 

Due  to  possible  difficulties  in  properly  removing  the  regional  and  local  geologically 
derived  fields  from  a  magnetic  survey,  we  also  allow  for  a  dc-shift  in  our  dipole  model; 
this  means  that  seven  parameters  are  required  to  characterize  each  anomaly. 


Inversion  goal:  Find  the  location  x  =  ( x,y,z ),  dipole  moment  m  =  (mX)  my.  mz) 
and  dc-shift  d,  that  produces  the  best  fit  to  a  set  of  TV  total-field  observations 
{do6s(xt),f  =  1,...,TV}. _ 
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Figure  5:  Spheroid  dimensions  that  can  produce  the  same  dipole  as  a  105  mm  shell  at  45° 
inclination  to  the  Earth’s  field.  The  dimensions  and  orientation  angles  corresponding  to 
the  labelled  points  are  shown  in  Table  1. 


We  follow  the  usual  convention  in  inverse  problems  and  use  m  to  specify  the  seven 
model  parameters  (it  should  be  clear  from  the  context  whether  m  refers  to  the  model 
parameters  or  the  dipole  moment).  We  formulate  the  inverse  problem  as  a  weighted, 
least-squares  optimization  procedure  where  the  goal  is  to  solve 


The  objective  function  is 


min  (b(m) 
meftn 


(12) 


4>( m)  =  l(F(m)  -  do6s)TW(F(m)  -  do6s)  (13) 

where  F(m)  is  the  forward  model  that  produces  the  predicted  data,  do6s  is  the  observed 
data  and  W  is  the  data-weighting  matrix  (ideally  the  data  weighting  matrix  is  the  inverse 
covariance  matrix  of  the  errors).  We  assume  independent,  gaussian  errors,  so  that  the 
data-weighting  matrix  is  diagonal  with  entries  Wa  =  \J af .  The  objective  function  may 
then  be  rewritten  as 


N 


<Km)  =  5  E 


t= 1  L 


Fi(  m)-df 


obs 


1  2 


N 


=  ^EriH2 


(14) 


i= 1 
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m  ♦ 


where  rj(m)  is  the  i-th  normalized  residual,  and  we  are  assuming  that  there  are  N 
observations. 

3.1  Inversion  method 

We  use  the  interior-reflective  Newton  method  (Branch  et  al.,  1999)  to  find  the  best 
fitting  model  parameters.  The  algorithm  is  a  trust  region  method  which  means  that  the 
vector  Sfc  that  updates  the  current  model,  i.e., 

mfc+i  =mk+sk  (15) 

is  chosen  from  a  Quadratic  model  of  the  objective  function  within  a  region  where  this 
model  can  be  trusted  (hence  the  term  trust  region  method).  A  quadratic  approximation 
to  <^>(m)  is 


<f>(mk  +  sfc)  «  <f>(mk)  +  V(f>(mk)sk  +  ^sJ’vV(mfc)sfc 
For  a  least  squares  problem  we  have 


(16) 


V0(m)  =  JT(m)  r(m)  (17) 

where  J  is  the  Jacobian  matrix  which  expresses  the  gradient  of  the  residuals  with  respect 
to  each  model  parameter, 


and  the  Hessian  matrix,  H ,  is  given  by 


drj{  m) 
drrij 


(18) 


N 

ff(m)  =  V2(f>(m)  =  JT(m)7(m)  +^ri(m)V2rj(m)  (19) 

Z=1 

The  Hessian  can  be  difficult  to  calculate  because  it  requires  second-order  derivatives. 
In  least-squares  problems  a  Gauss-Newton  approximation  is  often  used,  whereby  the 
second  term  in  the  above  equation  is  ignored.  This  is  a  good  approximation  when  the 
residuals  are  small  and  means  that  the  Hessian  can  be  obtained  essentially  for  free  once 
the  Jacobian  has  been  calculated. 

The  trust  region  method  involves  finding  the  s  (we  drop  the  subscripts  denoting 
iteration  number  and  the  dependence  on  m  for  clarity),  that  satisfies 

min  jisTJTJs  +  sTJTrj  such  that  ||Ds||  <  A  (20) 

where  D  is  a  diagonal  scaling  matrix  and  A  is  a  positive  scalar  that  controls  the  trust- 
region  size.  Without  the  trust  region  constraint  the  solution  of  the  above  equation  can 
be  found  from  the  normal  equations 


JTJs  =  -JT  r  (21) 

The  usual  trust-region  approach  is  the  Levenburg-Marquart  algorithm  where  we  solve 

(JTJ  +  A/)s  =  -JT  r  (22) 
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The  positive  parameter  A  ensures  that  the  solution  meets  the  constraint  ||Ds||  <  A  and 

also  improves  the  condition  number  of  the  normal  equations. 

Instead  of  minimizing  Equation  (20)  over  the  full  space,  the  algorithm  we  use  min¬ 
imizes  over  a  two-dimensional  subspace,  span(si,s2).  The  first  basis  of  the  subspace  is 
chosen  as  the  steepest  descent  direction 

si  =  -  JTr  (23) 

and  the  second  is  obtained  by  applying  preconditioned  conjugate  gradients  (PCG)  to 
Equation  (21).  The  PCG  iterations  are  terminated  when  a  negative  curvature  direction 
is  discovered,  i.e., 


s^JTJs2  <  0  (24) 

or  until  the  residual  is  sufficiently  small;  this  would  be  an  approximate  Gauss-Newton 
step.  The  philosophy  behind  the  method  is  to  force  global  convergence  (by  the  steep¬ 
est  descent  or  negative  curvature)  while  still  achieving  rapid  local  convergence  (by  the 
Gauss-Newton  step). 

3.1.1  Bounds  on  model  parameters 

The  Interior-reflective  Newton  Method  allows  bounds  to  be  placed  on  the  model  param¬ 
eters,  if  required.  We  impose  bounds  on  the  spatial  location  of  the  dipole, 

—2<x<2,  -2  <y  <2  and  -  4  <  z  <  0  (25) 

but  do  not  bound  the  dipole  moment  or  dc-shift.  With  lower  and  upper  bounds  of  [/, ,  u,] , 
the  minimization  problem  of  Equation  (12)  is  changed  to 

:  li  ^  mi  ^  M  (26) 

This  constrained  minimization  problem  can  be  solved  using  a  modification  of  the  for¬ 
mulation  presented  above  (Branch  et  al.,  1999). 

3.1.2  Starting  model 

Due  to  the  possible  presence  of  local  minima  it  is  important  to  obtain  a  good  first  guess 
of  the  model  parameters.  To  achieve  this  goal,  we  use  the  approach  of  McFee  and  Das 
(1986),  whereby  the  dipole  parameters  are  analytically  obtained  from  the  z-component 
of  the  vector  field.  Typically,  the  data  available  comprise  the  total-field  and  we  first  need 
to  estimate  the  z-component.  Most  anomalies  from  ordnance  have  peak  amplitudes  of 
less  than  1,000  nT,  whereas  the  Earth’s  field  is  on  the  order  of  50,000  nT.  The  anomalous 
total  field  ABt  is  then  approximately, 

A Bt  ~  axABx  +  ayABy  +  azABz  (27) 

where  a  =  ( ax,ay,az )  is  a  unit  vector  in  the  direction  of  the  Earth’s  field  and  AB  = 
(A Bx,  ABy,  ABZ)  is  the  vector  magnetic  anomaly.  Almost  all  of  the  continental  USA 
and  Canada  has  inclinations  exceeding  60°  so  that  az  >  0.85  and  to  a  first  approximation 


In  regions  where  the  z-component  of  the  field  is  not  so  dominant,  A Bz  can  be  obtained 
from  ABt  using  well  known  Fourier  transformation  techniques  (Blakely,  1996). 

The  McFee  and  Das  (1986)  technique  involves  first  locating  the  maximum  and  min¬ 
imum  values  of  the  field.  Then; 

•  The  line  joining  these  two  points  gives  the  azimuth  of  the  dipole; 

•  The  distance  between  the  two  points  is  related  to  the  depth  of  the  dipole; 

•  The  locations  of  the  extremum  and  their  relative  amplitude  can  be  used  to  infer 
the  horizontal  location  (assumed  to  lie  on  the  line  joining  the  points); 

•  The  relative  amplitudes  of  the  extremum  are  related  to  the  dipole  dip  angle; 

•  The  peak  anomaly  amplitude  is  proportional  to  the  dipole  magnitude;  and 

•  Comparing  the  fit  of  the  dipole  to  the  observed  data,  estimates  the  dc-shift. 

For  a  nearly  vertical  dipole  there  is  only  a  well  defined  maximum  (or  minimum). 
In  that  case,  the  maximum  gives  the  horizontal  location,  the  half-width  relates  to  the 
depth  and  the  maximum  amplitude  scales  with  the  magnitude  of  the  vertical  dipole. 


3.1.3  Data  weighting  matrix 

For  the  least-squares  data  fitting,  we  weight  each  data-point  dfs,  with  an  estimate  of 
the  standard  deviation  at  that  point,  crl,  through  the  data-weighting  matrix, 


(e  +  7Ko6s|)2 


(29) 


The  parameter  7  is  a  percentage  error  that  down-weights  the  influence  of  the  large  data 
values;  typically  we  chose  1%  <  7  <  5%.  The  parameter  e  ensures  that  the  denominator 
in  the  above  equation  does  not  get  too  close  to  zero  when  d°^s  approaches  zero;  typically 
we  choose  1  nT. 


3.1.4  Model  parameter  scaling 

If  different  model  parameters  have  significantly  different  scale  lengths  and  typical  values 
(e.g.  position  1  m;  dipole  moment  0.01  Am2),  then  small  favorable  changes  in  some 
model  parameters  may  be  hidden  by  the  large  changes  required  in  others.  The  ideal  is  to 
have  a  non-dimensionalized  model  space,  where  unit  changes  in  each  model  parameter 
have  about  the  same  impact  on  the  objective  function.  We  therefore  scale  each  model 
parameter  by  dividing  by  the  following  values, 

xtyp  =  (1, 1,  l)m,  mtyp  =  (0.1, 0.1, 0.1) Am.2  and  dtyp  =  5n T  (30) 

4  APPLICATION  OF  THE  INVERSION  METHOD 

Currently,  the  inversion  algorithms  are  implemented  in  Matlab  and  comprise  part  of  a 
GUI  driven  toolbox  for  geophysical  processing  developed  at  The  University  of  British 
Columbia.  Anomalies  are  first  identified  in  a  data-set  using  a  combination  of  automatic 
and  manual  interpretation.  Then,  the  inversion  algorithm  is  applied  to  a  user  specified 
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X  =  33.43  (±0.010)  m 
Y  =  10.71  (±0.010)  m 
Z  =  -0.25  (±0.008)  m 
Moment  =  0.0514  (±0.0017)  Am2 
Azimuth  =  -98.6°  (±6.1°) 

Dip  =73.2°  (±1.7°) 

CorrCoeff  =  0.970 
Angle  =  152.8°  (±1.7°) 

Successful  Fit 


Figure  6:  Example  result  of  dipole  moment  inversion. 


area  about  each  anomaly  (typically  around  2  x  2  m2  to  3  x  3  m2).  The  algorithm  is 
terminated  when  the  objective  function  decreases  by  less  than  a  specified  tolerance, 
or  when  the  gradient  falls  below  a  specified  value.  Both  these  parameters  are  user 
controlled. 

Summary  reports  can  be  generated,  and  the  results  visually  inspected  as  demon¬ 
strated  in  Figure  6.  We  will  illustrate  the  performance  of  the  inversion  algorithm  on  a 
large  magnetics  survey  collected  as  part  of  an  clearance  project  in  the  Helana  Valley, 
Montana. 

4.1  The  Guthrie  Road,  Montana  dataset 

Once  open  fields,  portions  of  the  Helena  Valley  were  used  in  the  1950s  for  military 
training  by  Montana  Army  National  Guard  (MTARNG).  The  National  Guard  Bureau 
funded  the  remediation  project  which  was  divided  into  two  phases  focusing  on  different 
impact  areas:  Diamond  Springs  and  Guthrie  Road.  We  consider  the  data  collected  in 
the  Guthrie  Road  area  in  the  summer  of  1998.  Geophysical  Technology  Limited  (GTL) 
of  Armidale,  Australia,  conducted  a  magnetometer  survey  of  the  area  using  an  all-terrain 
vehicle  towing  an  array  of  eight  cesium  vapor  magnetometers  (Clark  et  ah,  1999).  The 
system  was  equipped  with  a  real-time  global  positioning  system  that  allowed  anomalies 
to  be  positioned  to  within  about  20  centimeters. 

GTL  found  840  anomalies  that  were  tagged  as  potential  UXO,  and  these  were  exca¬ 
vated  in  the  summer  of  1999.  We  analyze  822  of  the  anomalies  (we  couldn’t  locate  the 
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Figure  7:  Some  of  the  81  mm  mortars  recovered  at  Guthrie  Rd. 


data  for  the  remaining  18  anomalies).  Table  2  partitions  the  dig-sheet  identifications 
into  six  categories;  76  mm  ordnance,  81  mm  mortars,  large  pieces  of  ordnance,  shrapnel, 
metallic  debris  and  geology  (including  anomalies  where  no  metallic  items  were  found). 

A  selection  of  the  81  mm  mortars  recovered  at  Guthrie  Rd  are  shown  in  Figure  7. 
Most  are  intact  and  have  not  been  substantially  deformed  by  the  force  of  impact.  How¬ 
ever,  there  are  a  number  of  deformed  but  intact  items.  These  tend  to  have  been  practice 
rounds  that  were  deformed  when  the  black  spotter  charge  detonated  on  impact.  There 
axe  also  a  number  of  items  that  have  been  split  into  pieces,  these  were  not  classified  as 

uxo. 

We  fitted  dipole  moments  using  a  percentage  error  of  7  =  5%  and  a  base  error  of 
e  =  1  n T.  The  entire  run  of  inversions  took  around  half  and  hour  on  a  fairly  standard 
desktop  computer  (Pentium  Pro  Processor  with  800  MHz  processor  and  1  Gb  RAM). 
We  use  the  correlation  coefficient  c  between  fitted  and  observed  data  to  determine  the 
goodness  of  fit.  We  found  that  85,  or  around  10%  of  anomalies  had  fits  with  c  <  0.7  and 
162,  or  around  20%  had  c  <  0.8.  A  break-down  according  to  the  dig-sheets  of  the  failed 
fits  is  given  in  Table  2.  Most  of  the  failed  fits  are  from  items  identified  as  shrapnel, 
geology  or  metallic  debris,  with  only  one  ordnance  item  having  a  failed  fit.  We  did  not 
visually  inspect  every  failure,  but  those  we  did  inspect  generally  had  some  type  of  data 
issue.  These  included  overlapping  acquisition  lines  that  we  couldn’t  remove  without 
undue  effort;  levelling  problems  between  adjacent  lines;  and  overlapping  anomalies.  In 
no  case  did  the  inversion  algorithm  fail  where  the  data  quality  was  good.  In  the  next 
sections  we  will  describe  methods  for  discrimination  and  identification  and  apply  them 
to  this  data-set. 


5  UXO  DISCRIMINATION 


The  last  section  described  how  dipole  moments  and  position  are  fit  to  each  anomaly 
by  our  inversion  routine.  If  the  data-fit  is  unsatisfactory  the  anomaly  is  marked  as 
such  and  in  a  production  clean-up  would  be  excavated  as  potential  UXO.  From  our 
viewpoint  as  data  interpreters  we  are  unable  to  provide  any  reliable  information  on  the 
target  characteristics.  On  the  other  hand,  if  the  fit  is  satisfactory  we  attempt  to  use 
the  recovered  dipole  moment  to  discriminate  UXO’s  from  non-UXO’s.  For  those  items 
classified  as  potential  UXO  the  additional  step  of  ordnance  identification  is  attempted 
as  described  in  the  following  section. 

Our  method  for  ordnance  discrimination  is  based  on  that  described  earlier  in  this 
report,  and  is  due  to  Nelson  et  al.  (1998).  Briefly,  the  method  relies  on  the  following 
two  factors; 


1.  Shock  demagnetization  causes  the  remnant  magnetism  of  an  ordnance  to  be  vir¬ 
tually  eliminated; 

2.  The  induced  magnetization  in  a  spheroid  with  e  <  7  is  constrained  to  lie  within 
about  60°  of  the  Earth’s  field  (Figure  3). 


With  a  recovered  dipole  moment  of  m,  our  discrimination  routine  is  to  calculate  the 
angle  between  it  and  the  Earth’s  field,  B0, 


ip  =  arccos( 


m  •  B0 

iHiiiBoir 


(31) 


We  can  never  be  certain  that  all  remnant  magnetization  has  been  erased  by  the  shock 
of  impact.  Additionally,  non-ordnance  items,  whether  remnantly  magnetized  or  not, 
may  have  dipole  angles  close  to  the  Earth’s  field.  Consequently,  it  is  not  possible  to 
make  a  definitive  ordnance/non-ordnance  distinction  for  all  anomalies.  Our  preferred 
discrimination  method  is  to  rank  items  according  to  the  likelihood  they  are  due  to 
ordnance.  We  do  this  on  the  basis  of  the  angle  ip;  the  items  with  smaller  ip  are  ranked 
higher  up  the  list.  In  a  clean-up  operation  one  then  excavates  according  to  the  list.  If 
resources  are  available  to  dig  50%  of  the  anomalies  then  the  first  50%  of  items  in  the 
list  are  excavated.  As  a  quality  control  measure  a  certain  number  of  the  lower  ranked 
items  should  also  be  excavated. 


5.1  The  Guthrie  Road,  Montana  dataset 

In  Figure  (8a)  we  plot  the  dipole  magnitude  and  angle  relative  to  the  Earth’s  field  of  all 
the  non-ordnance  anomalies.  Most  of  the  dipoles  lie  within  ±90°  of  the  Earth’s  field  but 
a  large  number  have  angles  exceeding  60°.  On  the  other  hand,  the  dipole  angles  for  the 
ordnance  and  large  ordnance  fragments  almost  all  lie  within  ±60°  of  the  Earth’s  field 
Figure  (8b).  These  results  indicate  that  significant  shock  demagnetization  has  occurred 
to  the  ordnance.  There  is  only  one  significant  outlier,  with  an  angle  of  around  70°. 

In  Figure  (9a)  we  give  a  normalized  histogram  of  the  angles  for  the  different  categories 
of  items.  There  is  a  clear  concentration  of  ordnance  with  dipole  angles  close  to  the 
Earth’s  field  (we  include  the  significant  sized  ordnance  fragments  in  this  category). 
Shrapnel  tends  to  have  relatively  low  angles,  but  there  are  a  significant  number  of  items 
with  larger  angles.  There  are  two  simple  mechanisms  that  could  explain  this  distribution: 
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Figure  8:  (a)  Polar  plot  of  dipole  angle  and  magnitude  (Am2)  of  non-ordnance  items; 
(b)  Same  for  ordnance  (cross)  and  significant  ordnance  fragments  (dots). 


Figure  9.  Breakdown  of  dig-sheets  identifications  and  their  (a)  dipole  angle  relative  to 
the  Earth’s  magnetic  field;  and  (b)  dipole  magntitude. 


1.  Shock  demagnetization  sometimes  occurs  with  shrapnel  but  not  always; 

2.  Shrapnel  acquires  a  remnant  magnetization  after  impact,  in  the  presence  of  the 
Earth’s  field. 

Metallic  debris  and  geological  anomalies  can  have  significant  remnant  magnetism, 
with  a  definite  peak  in  the  distribution  between  60°  and  100°.  We  can  think  of  no 
plausible  explanation  for  this  peak  in  the  distribution  of  dipole  angles. 

The  results  of  discrimination  using  different  cut-off  angles,  for  intact  ordnance  items 
are  summarized  in  Table  3,  and  for  ordnance  items  and  large  fragments  in  Table  4.  The 
tables  give  the  number  of  items  to  dig  (including  the  85  anomalies  with  poor  fits),  the 
number  left  in  the  ground,  the  ordnance  found,  the  ordnance  left  and  the  number  of 
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Figure  10:  Yield  of  ordnance  with  number  of  holes  dug  for  the  three  discrimination 
methods.  For  all  methods  it  is  assumed  the  85  anomalies  with  poor  fits  are  excavated 
first. 

false  alarms.  Figure  (10)  plots  the  proportion  of  items  recovered  versus  the  number  of 
holes  that  need  to  be  excavated.  Once  the  85  anomalies  with  poor  fits  have  been  dug, 
there  is  a  high  yield  of  ordnance  as  the  number  of  holes  increases.  However,  once  about 
80%  of  items  have  been  recovered  there  is  a  noticeable  slow-down  in  the  ordnance  yield; 
a  lot  of  holes  need  to  be  dug  to  recover  the  last  few  ordnance  items.  A  total  of  560 
(68%)  of  holes  are  required  to  recover  all  ordnance  items;  this  number  also  recovers  all 
the  large  ordnance  pieces. 

Many  of  the  dipole  fits  for  shrapnel  and  metallic  debris  have  moments  with  low 
magnitude  (Figure  9b).  None  of  the  intact  ordnance  items  has  a  dipole  magnitude  less 
than  0.05  Am2.  However,  seven  of  the  large  ordnance  fragments  have  dipole  moments 
less  than  this  value.  The  results  of  rejecting  items  as  non-ordnance  based  on  the  dipole 
magnitude  (as  well  as  a  cutoff  angle  of  75°)  are  summarized  in  Tables  5  (for  ordnance 
only)  and  6  (for  large  fragments  as  well).  With  a  dipole  cutoff  of  0.05  Am2  no  intact 
ordnance  items  are  missed  and  we  can  leave  379  items  in  the  ground  (46%  saving).  If  we 
make  the  reasonable  assumption  (for  this  area)  that  the  smallest  item  likely  to  be  found 
is  a  60  mm  mortar,  the  0.05  Am2  cutoff  is  quite  sensible;  the  minimum  anomaly  from 
a  60  mm  mortar  is  around  0.055  Am2.  Ranking  items  by  their  magnitude  results  in  a 
lower  initial  yield  of  ordnance  (Figure  10)  than  ranking  them  by  angle.  However,  the 
magnitude  ranking  does  not  suffer  as  significant  a  slow-down  in  yield,  with  a  crossover 
occurring  once  about  80%  of  items  have  been  recovered.  A  total  of  443  (54%)  of  holes 
need  to  be  excavated  to  recover  all  the  ordnance. 
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We  try  one  additional  discrimination  method.  It  uses  the  idea  of  remnant  magnetism 
and  requires  that  we  have  some  information  on  the  ordnance  items  likely  to  be  found 
in  the  area.  For  the  Guthrie  Rd  data  we  assume  that  60  and  81  mm  mortars  and  76, 
90,  105  and  155  mm  projectiles  may  be  present  in  the  area.  Our  method  is  described 
in  more  detail  in  the  next  section.  Briefly,  we  work  out  the  minimum  percentage  of 
remnant  magnetism  that  is  required  to  make  the  recovered  dipole  moment  match  of 
the  ordnance  items.  We  then  rank  the  items  according  to  the  percentage  of  remnant 
magnetism,  with  lower  percentages  higher  up  the  list.  The  method  does  not  apply  to 
large  ordnance  fragments  as  these  are  not  considered  in  our  modelling. 

The  results  of  discriminating  by  remnant  magnetism  are  summarized  in  Table  7. 
Note,  the  we  also  incorporate  the  angle  constraint  (<  75°)  and  the  magnitude  con¬ 
straint  (>  0.05  Am2).  Allowing  up  to  half  of  the  dipole  anomaly  to  be  due  to  remnant 
magnetization,  we  can  identify  all  ordnance  items  and  leave  407  anomalies  in  the  ground 
(49.5  %  saving).  The  ordnance  versus  number  of  holes  curve  for  remnant  magnetization 
((Figure  10)  reveals  a  high  yield  of  ordnance  with  increasing  number  of  excavations. 
There  is  only  a  very  slight  slowdown  in  yield  once  about  95%  of  ordnance  items  have 
been  found.  Once  402  (49%)  of  anomalies  are  excavated  all  ordnance  items  are  recovered. 

6  ORDNANCE  IDENTIFICATION 

For  those  anomalies  classified  as  potential  ordnance,  the  next  step  is  to  try  to  identify 
the  ordnance  type. 

6.1  IDENTIFICATION  USING  THE  RECOVERED  DIPOLE  MOMENT 

We  established  in  Appendix  B  and  an  earlier  part  of  the  report  that  a  given  dipole 
anomaly  could  have  been  produced  by  an  infinite  number  of  spheroids  (Figure  5).  How¬ 
ever,  for  any  given  site,  there  are  only  a  finite  number  of  ordnance  types  that  are  likely 
to  be  found.  This  is  the  key  observation  that  allows  us  to  identify  ordnance  type  in  the 
presence  of  ambiguity. 

The  dipole  induced  in  a  spheroid  depends  on  the  angle  between  the  spheroids  semi¬ 
major  axis  and  the  Earth’s  field.  There  is  no  azimuthal  dependence  due  to  symmetry. 
As  the  ordnance  is  rotated  about  the  Earth’s  field  the  induced  dipole  moment  will  trace 
out  a  curve  such  as  that  shown  in  Figure  (4).  This  represents  all  dipole  moments  that 
can  be  produced  by  the  given  ordnance  item  in  the  absence  of  remnant  magnetism. 
This  mapping  procedure  can  be  repeated  for  all  ordnance  items  suspected  of  occurring 
in  the  area.  With  N0  ordnance  items,  a  series  of  curves,  {mk(0),  k  =  1, . . .  ,iV0},  will 
be  generated.  For  a  particular  recovered  dipole  moment  m,  we  calculate  the  minimum 
distance  between  it  and  each  of  the  dipole  curves 

Am*  =  min  llm  -  im(0)||  (32) 

O<0<2tt  1  wn  V  ' 

We  assume  that  any  discrepancy  is  due  to  remnant  magnetism  which,  as  a  percentage 
of  the  observed  moment,  will  be 
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Figure  11:  Dipole  moments  of  the  Guthrie  Road  ordnance  items  for  (a)  76  mm  projec¬ 
tiles;  and  (b)  81  mm  mortars. 


This  allows  us  to  rank  the  ordnance  items  according  to  the  likelihood  they  produced 
the  observed  anomaly.  The  most  likely  item  is  that  with  the  lowest  percentage  remnant 
magnetism. 

The  success  of  the  procedure  will  depend  on; 

1.  The  accuracy  of  the  recovered  dipole  moment  (see  the  error  analysis  section); 

2.  The  orientation  of  the  ordnance,  because  some  orientations  produce  dipoles  that 
are  close  to  dipoles  potentially  produced  by  other  ordnance  items; 

3.  How  closely  the  spheroid  approximation  matches  that  for  real  ordnance;  and 

4.  How  much  remnant  magnetism  the  ordnance  maintains  after  the  shock  of  impact; 

6.1.1  The  Guthrie  Road,  Montana  dataset 

The  dipole  moments  for  those  anomalies  due  to  76  mm  projectiles  are  shown  in  Fig¬ 
ure  (11a)  and  for  81  mm  mortars  in  Figure  (lib).  Misidentifications  for  both  types  of 
ordnance  appear  likely  due  to  the  overlapping  curves  for  the  difference  ordnance  types. 
In  particular  the  moments  for  almost  all  the  76  mm  projectiles  are  in  a  region  where  the 
60,  76  and  81  mm  calibre  ordnances  have  similar  dipole  moments.  For  the  81  mm  mor¬ 
tars  about  half  the  moments  are  in  a  region  with  only  the  curve  for  the  81  mm  nearby, 
while  the  other  half  are  in  the  more  complex  part  of  the  plot.  Our  dipole  identification 
routine  involves  finding  the  orndance  curve  that  is  closest  to  each  recovered  moment. 

Dipole  identification  results  for  76  mm  projectiles  are  given  in  (Figure  12a)  and  for 
81  mm  mortars  in  Figure  (12b).  In  19  out  of  34  cases  (56%)  the  76  mm  projectile  is 
identified  correctly.  The  most  common  case  of  misidentification  is  as  a  60  mm  mortar; 
from  Figure  (4)  we  see  that  this  is  due  to  overlapping  moments  between  the  60  and  76 
mm  ordnances.  For  81  mm  mortars  the  identification  routine  is  correct  in  25  of  48  cases 
(52%),  so  that  overall  we  achieve  correct  identifications  54%  of  the  time.  The  most  com¬ 
mon  misidentification  for  the  81  mm  mortar  is  as  a  76  mm  projectile  (13  times).  Some 
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Figure  12:  Identification  results  for  the  Guthrie  Road  data  using  both  the  dipole  and 
dipole  plus  octupole  methods  for  (a)  76  mm  projectiles;  and  (b)  81  mm  mortars. 


of  the  errors  in  misidentification  could  relate  to  the  deformation  of  certain  of  the  81 
mm  practice  mortars  that  occurred  when  the  spotter  charge  detonated  (Figure  13).  A 
spheroid  is  unlikely  to  give  a  very  good  approximation  to  this  item.  Unfortunately,  while 
recovered  ordnance  items  were  retained  after  excavation  no  records  of  the  anomaly  num¬ 
ber  were  kept.  Thus,  we  cannot  determine  which  recovered  dipole  moments  correspond 
to  deformed  ordnance. 

Assuming  we  are  using  the  remnant  magnetization  method  for  discrimination  and  use 
a  cutoff  of  50%,  there  will  be  235  false  alarms;  items  incorrectly  identified  as  ordnance. 
Applying  the  identification  method  to  these  false  alarms  (Figure  14a)  is  most  likely  to 
cause  a  misidentification  as  a  60  mm  mortar  (the  moments  are  generally  quite  low). 
Metallic  debris  is  more  likely  to  give  a  false  identification  as  a  large  ordnance  item  than 
any  of  the  other  categories  of  items  (shrapnel,  large  ordnance  pieces  and  geology). 

6.1.2  Seeded  test  site  at  The  Former  Fort  Ord,  California 

The  former  Fort  Ord  is  located  near  Monterey  Bay  in  northwestern  Monterey  County, 
California.  Since  1917,  portions  of  Fort  Ord  were  used  by  Army  units  for  maneuvers, 
target  ranges,  and  other  purposes.  Investigation  and  ordnance  removal  actions  have 
been  conducted  at  the  former  Fort  Ord  since  1994  when  the  base  closed.  In  November 
1998,  the  United  States  Department  of  the  Army  instigated  an  Ordnance  Detection  and 
Discrimination  Study  (ODDS)  at  Fort  Ord  (Asch  and  Staes,  2001). 

The  purpose  of  the  Seeded  Test  was  to  evaluate  the  detection  and  discrimination 
capabilities  of  selected  geophysical  instruments  when  ordnance  items  are  buried  or 
’’seeded”  below  the  ground  surface  in  a  local  site-specific  geologic  conditions.  The  seeded 
test  area  consists  of  five  separate,  adjacent  grids,  which  were  cleared  of  anomalies  to  pro¬ 
vide  a  controlled  testing  environment.  Two  grids  are  ’’known”  plots  where  the  items, 
depths  and  orientations  of  items  are  made  known  to  all  participants.  It  is  these  two 
known  sites  (K1  and  K2)  that  we  analyze  here. 

Seventeen  different  types  of  ordnance  items  were  buried  in  the  seeded  test  plots. 
Some  ordnance  scrap  and  other  items  were  also  seeded  to  allow  evaluation  of  discrimi- 
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Figure  13:  One  of  the  practice  81  mm  mortars  deformed  due  to  the  detonation  of  the 
spotting  charge. 


nation  capabilities.  Additional  areas  of  the  plots  were  excavated  and  back-filled  without 
any  buried  items. 

For  the  K1  5  of  the  7  ’’empty  hole”  items  (618,  619,  621  to  623)  did  not  produce 
acceptable  fits,  so  they  would  not  have  been  declared  as  potential  ordnance.  The  other 
two  ’’empty  holes”  had  poorly  fitting  dipoles  and  would  have  been  identified  as  scrap. 
A  fragment  (595)  and  a  signal  illumination  flare  (612)  also  produced  unacceptable  fits. 
Discrimination  and  identification  results  for  the  remaining  items  are  summarized  in 
Table  9.  Two  of  the  four  fragments  were  correctly  identified  as  non-ordnance  (using 
a  discrimination  cutoff  of  60°).  Six  of  the  twelve  35  mm  projectiles  were  incorrectly 
classified  as  scrap;  this  may  be  due  to  remnant  magnetism.  Five  of  the  remaining  35 
mm  projectiles  had  the  correct  identification  ranked  in  the  top  three.  One  of  the  2.36” 
rockets  was  classified  as  scrap,  the  other  three  had  a  2.36”  identification  ranked  in  the 
top  three  (one  was  correct).  The  six  signal  illumination  flares  could  not  be  identified 
correctly,  neither  could  any  of  the  grenades.  Overall,  the  discrimination /identification 
results  for  K1  are  not  very  impressive.  We  believe  part  of  the  reason  could  be  due  to 
remnant  magnetism  (the  items  were  not  shock  demagnetized  prior  to  emplacement).  It 
may  also  be  a  reflection  of  the  generally  poor  performance  of  magnetics  for  small  calibre 
items  (Klaff  et  al.,  2002). 

For  the  K2  plot,  good  dipole  fits  could  not  be  obtained  for  four  items,  555,  556  (both, 
25  mm  projectiles  at  46  cm),  569  (fragment  at  30  cm)  and  574  (signals  at  30  cm).  The 
results  for  the  rest  of  the  items  are  summarized  in  Table  10.  There  were  seventeen  items 
correctly  identified  as  ordnance.  Six  of  these  were  identified  as  the  correct  ordnance 
type  and  five  others  had  the  correct  ordnance  listed  in  the  top  three  possibilities.  The 
other  six  identifications  were  a  long  way  off  and  we  believe  that  several  of  these  were 
remnantly  magnetized;  in  particular  all  the  stokes  mortars  and  the  hand  grenade. 
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Figure  14:  Identification  results  for  the  Guthrie  Road  data  for  (a)  dipole  method;  and 
(b)  dipole  plus  octupole  method. 

6.2  IDENTIFICATION  USING  THE  DIPOLE  AND  OCTUPOLE 

Magnetic  anomalies  from  ordnance  are  usually  dominated  by  the  dipole  component; 
however,  with  high  quality  data  we  may  be  able  to  obtain  useable  information  from 
the  octupole  term.  Due  to  the  dipole  dominance  it  would  probably  not  be  a  good  idea 
to  attempt  a  direct  inversion  for  ordnance  dimensions  (local  minima,  non-uniqueness 
etc).  Rather,  a  better  method  is  to  obtain  an  optimum  fit  for  each  ordnance  item  in 
our  library,  and  declare  the  ordnance  with  the  best  fit  as  the  most  likely  source  of  the 
anomaly.  For  each  item  in  the  library  this  procedure  involves  determination  of  six  model 
parameters  ( x ,  y ,  z,  6 , <£,  d),  where  ( x ,  y,  z)  is  location,  6  is  the  ordnance  dip  angle  relative 
to  horizontal,  <t>  is  the  azimuth  angle  relative  to  geographic  North  and  d  is  the  dc-shift. 

We  again  use  a  least-squares  formulation  and  the  interior-reflective  Newton  method. 
The  data-weighting  and  parameter  bounds  are  identical  to  that  used  for  the  dipole 
inversion  (i.e.  we  only  bound  the  spatial  location).  For  scale  factors  we  use 

xtyp  =  (1,1,  l)m,  ( Qtyp ,  <f>tyP)  =  (?r/2,  7t/2)  radians  and  dtyp  =  5n T  (34) 

To  determine  the  starting  model,  we  first  do  a  dipole  inversion  and  go  through  the  dipole 
identification  process.  For  each  ordnance  item,  the  identification  process  will  return  the 
dip  angle  6  and  azimuth  0  that  achieved  the  closest  fit  to  the  dipole.  These  can  be 
used  as  the  starting  orientation  with  the  dipole  location  used  as  the  starting  location, 
similarly  for  the  dc-shift.  Rather  than  doing  an  inversion  for  every  item  in  the  library 
we  only  consider  those  3-5  ordnance  items  with  the  best  dipole  fits. 

6.2.1  The  Guthrie  Road ,  Montana  dataset 

We  only  apply  the  dipole/octupole  identification  procedure  to  those  anomalies  flagged  as 
potential  UXO  by  the  remnant  magnetization  technique  (with  a  cutoff  of  50%  remnant 
magnetism).  The  quality  of  the  fits  obtained  for  the  optimum  ordnance  were  around 
the  same  as  the  dipole  fits;  sometimes  they  were  slightly  worse  and  sometimes  slightly 
better. 
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Identification  results  for  the  dipole/octupole  method  for  76  mm  projectiles  are  given 
in  Figure  (12a)  and  for  81  mm  mortars  in  Figure  (12b).  For  76  mm  projectiles  the 
correct  identification  is  achieved  for  20  of  34  cases,  one  better  than  the  dipole  method. 
However,  for  the  81  mm  mortar  correct  identification  is  only  achieved  for  22  of  48  cases, 
three  less  than  the  dipole  method. 

Misidentification  results  for  false  alarms  are  given  in  Figure  (14b).  The  distribution 
of  items  is  very  similar  to  that  achieved  with  the  dipole  method. 

We  believe  that  there  are  three  reasons  why  the  dipole/octupole  identification  method 
does  not  improve  the  original  dipole  identification  method. 

1.  The  octupole  response  is  quite  small  and  does  not  provide  significant  additional 
constraints  on  the  the  geometry; 

2.  The  equivalent  spheroid  is  not  an  accurate  enough  approximation  of  an  ordnance 
item  for  the  octupole  contribution  to  have  a  positive  effect;  and 

3.  The  ordnance  items  are  slightly  remnantly  magnetized  which  changes  the  relation¬ 
ship  between  ordnance  orientation  and  the  direction  of  magnetization. 

6.3  IDENTIFICATION  USING  THE  DIPOLE  AND  QUADRUPOLE 

Thus  far,  the  identification  routines  we  have  developed  have  relied  solely  on  a  spheroid 
model  of  ordnance.  However,  this  assumption  is  overly  simplistic  since  most  ordnance 
items  do  not  have  that  degree  of  symmetry.  In  particular,  the  front  and  tail  ends  of 
ordnance  tend  to  be  different:  the  front  is  often  tapered  to  provide  better  aerodynamics 
and  target  penetration,  while  the  tail  end  is  often  flat  (the  tails  of  mortars  don’t  matter 
because  they  are  usually  made  of  aluminum).  Therefore,  a  more  realistic  model  of 
ordnance  is  an  axially  symmetric  body  without  front-back  symmetry.  Such  a  body  has 
a  multipole  expansion  with  both  a  dipole  term  and  a  non-zero  quadrupole  term  (recall 
that  the  quadrupole  was  zero  for  a  spheroid  due  to  symmetry). 

Our  intent  with  this  new  identification  method  is  to  try  to  use  the  information  in  the 
quadrupole  to  constrain  the  orientation  of  the  body  (</>,  0).  As  we  discovered  when  using 
the  dipole  method,  the  principal  source  of  ambiguity  in  identification  is  our  inability 
to  constrain  the  orientation  of  the  body.  In  particular,  when  the  object’s  dimension  is 
changed,  one  can  compensate  by  changing  the  object  orientation  and  produce  the  same 
or  a  very  similar  dipole  moment.  If  we  could  constrain  the  orientation  then  we  could 
better  constrain  the  ordnance  dimension. 

The  basis  of  the  idea  is  that  there  will  be  no  quadrupole  contributions  perpendicular 
to  the  axis  of  symmetry  but  there  will  be  contributions  along  the  symmetry  axis.  We 
therefore  construct  a  model  with  these  characteristics  and  use  an  inversion  procedure 
to  recover  a  model  with  the  best  orientation.  Concomitant  with  this  process  we  also 
recover  parameters  that  depend  on  the  shape  of  the  body,  as  we  now  explain. 

The  dipole  moment  of  an  axially  symmetric  body  can  be  written  in  the  form, 


m  =  Mo  1AtFAB0  (35) 

where  A  is  the  usual  Euler  rotation  tensor  and  F  is  a  3  x  3  diagonal  matrix  with 
(F2,F2,Fs)  along  the  diagonal.  We  are  assuming  axial-symmetry  so  that  F\  =  F2. 
The  demagnetization  factors,  Fit  given  in  the  above  equation,  are  VFi  times  those 
previously  given  in  Equation  2,  where  V  is  the  volume  of  the  body.  These  volume 
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modified  demagnetization  factors  contain  all  the  information  on  the  body’s  dimensions 
that  we  are  attempting  to  recover. 

We  are  assuming  the  body  is  axially  symmetric  (about  the  z-axis),  therefore  the 
quadrupole  tensor  components  in  the  x-y  plane  are  zero  (see  Appendix  A  for  a  deriva¬ 
tion).  The  body  is  asymmetric  about  the  x-y  plane  so  that  the  quadrupole  tensor  has 
non-zero  components  in  the  z-direction.  The  tensor  is  symmetric  Mq  =  and  is 
therefore  defined  by  the  equation, 


Mq  = 


0  0  mqi 

0  0  mq  2 

Tflq\  Tflq2  rflq^ 


(36) 


where  mq 2  and  mq 3  are  dependent  on  the  shape  and  magnetization  of  the  body. 

The  above  equation  is  in  body-centered  coordinates.  The  quadrupole  tensor  may  be 
rotated  into  the  Earth’s  frame  of  reference  using  the  following  equation, 


Mg  =  AMgAT  (37) 

The  magnetic  field  due  to  this  quadrupole  tensor  can  be  calculated  using  Equation  (A- 
11).  Ignoring  the  octupole  and  any  higher  order  moments,  full  specification  of  an  axially- 
symmetric  body  therefore  involves  the  following  11  parameters 


[x,  y,  z,  <j>,  e,  Ft,  F3,  mq\,mq2,  m,3)  d]  (38) 

where  we  have  again  allowed  a  dc-shift,  d. 

To  invert  for  the  11  parameter  model  we  use  the  same  least-squares  formulation  and 
Interior-reflective  Newton  method  described  previously.  The  parameter  scalings  are  the 
same  as  those  given  in  Equation  (34),  with  the  additional  parameters  scaled  as 

Ftyp  =  (0.01, 0.01)m3  and  mqtyp  =  (0.01,  .01,  .01)Am3  (39) 

Once  the  model  parameters  have  been  recovered  we  return  to  our  spheroid  model  of 
ordnance  and  find  the  object  dimensions  that  reproduce  the  recovered  volume  demag¬ 
netization  factors,  F2  and  F3.  This  is  straightforward  because  the  aspect  ratio,  e  solely 
determines  the  ratio  F2/F3,  while  the  diameter  and  aspect  ratio  together  determine  the 
magnitudes  of  F2  and  F3.  On  occasion  we  find  that  F2  converges  towards  zero  while  F3 
appears  to  converge  to  a  sensible  value.  In  that  case,  we  fix  the  aspect  ratio  at  a  typical 
value  3  <  e  <  5  and  just  solve  for  the  diameter. 


6.3.1  The  Guthrie  Road ,  Montana  dataset 

We  applied  the  quadrupole  inversion  method  to  all  737  anomalies  at  Guthrie  Rd  that 
had  good  dipole  fits.  Perhaps  self-evidently  the  fits  and  correlation  coefficients  for 
the  dipole/quadrupole  model  were  better  than  for  the  dipole  model  alone.  Figure  15 
summarizes  the  results  of  the  quadrupole  identification.  There  appears  to  be  some 
identification  ability  but  there  are  a  significant  number  of  ordnance  items  where  the 
recovered  dimensions  are  inaccurate.  In  particular,  there  is  a  tendency  to  either  recover 
too  low  an  aspect  ratio  and  hence  too  large  a  diameter  or  vice-versa. 

There  are  two  conclusions  that  can  be  drawn  from  the  quadrupole  discrimination 
results.  Firstly,  the  quadrupole  may  be  unable  to  provide  the  required  constraints  on 
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Figure  15:  Quadrupole  identification  results  at  Guthrie  Road;  (a)  Recovered  aspect- 
ratios  and  diameters;  and  (b)  Histogram  of  recovered  diameters. 


the  ordnance’s  orientation.  Alternatively,  there  may  be  enough  remnant  magnetization 
in  the  body  so  that  our  model  is  not  an  accurate  reflection  of  reality. 

With  these  issues  in  mind  we  trialed  one  last  method  for  identification  using  the 
quadrupole.  This  time  we  allow  the  dipole  moment  to  vary  independently  of  the 
quadrupole;  i.e.  in  place  of  F\  and  F2  in  Equation  (38)  we  solve  for  m  =  (7721,7712,7713). 
Thus  there  are  twelve  parameters  to  find.  Our  intent  is  to  determine  if  unmodelled 
quadrupole  components  bias  the  estimation  of  the  dipole  moment. 

We  applied  the  method  to  the  Guthrie  Road  data  and  then  used  the  recovered  dipole 
moment  to  achieve  ordnance  identification  in  the  usual  way.  We  found  that  the  method 
was  slightly  inferior  to  the  original  dipole  identification  routine.  In  particular  for  the 
81  mm  mortar  23  of  the  48  anomalies  were  identified  correctly  (compared  to  25  of  48), 
while  for  the  76  mm  projectile  only  16  of  34  cases  were  correct  (compared  to  20  of  34). 
Thus,  we  conclude  that  unmodelled  quadrupole  components  do  not  negatively  affect  our 
ability  to  identify  using  the  dipole  moment. 

7  ERROR  ANALYSIS 

The  dipole  identification  results  for  Guthrie  Road  indicate  that  around  55%  of  the  time 
we  correctly  identify  the  ordnance  type.  The  45%  misidentifications  must  be  due  to 
either  an  incorrect  model,  remnant  magnetism  or  uncertainty  in  the  recovered  dipole 
moments.  In  this  section  we  analyze  the  uncertainty  in  recovered  dipole  moment  using 
linearized  error  analysis. 

7.1  Estimating  the  covariance  matrix  of  the  model  parameters 

To  derive  the  covariance  matrix  of  the  model  parameters  we  follow  the  development 
given  in  Bard.  At  the  minimum  of  the  objective  function  the  gradient  is  zero 
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V0(m,d)  =  O  (40) 

where  V  is  with  respect  to  the  model  parameters.  If  we  change  the  data  to  d  +  dd  we 
will  change  the  objective  function  and  hence  the  position  of  the  minimum  will  vary  to 
m  +  dm.  The  gradient  of  the  new  objective  function  at  the  new  minimum  will  be  zero. 

V<?l>(m  +  dm,  d  +  dd)  =  0  (41) 

If  we  expand  this  gradient  in  a  Taylor  series  and  only  keep  first  order  terms  we  find 

r\ 

V0(m  +  dm,  d  +  dd)  =  V<£(  m,  d)  +  V2<p(m,  d)dm  +  — V<£(m,  d)dd  =  0  (42) 

Recognizing  that  the  V2  term  is  the  Hessian  and  rearranging  we  find 

dm  =  H"1^V^(m,d)dd  (43) 

The  covariance  matrix  of  the  model  parameters  is 

Vm  =  E[dmdmT]  (44) 

where  E[-]  is  the  expectation  operator.  Substituting  Equation  (43)  into  (44),  recognizing 
that  the  Hessian  and  gradient  terms  are  constants  (under  expectation)  and  rearranging 
terms  we  find 

Vm  =  H-1— V^(m,d)TE[ddddT]— V<^(m,d)H-1  (45) 

The  covariance  matrix  of  the  data  is  =  E[<Jd<5dT].  Assuming  the  data-weighting 
matrix  in  our  formulation  of  the  inverse  problem  is  the  same  as  the  inverse  covariance 
matrix  of  the  data,  the  Gauss-Newton  approximation  to  the  Hessian  is 


H  =  J  TV^J  (46) 

while 

^V0(m,d)  =  Vd-1J  (47) 

Substitution  of  Equations  (46)  and  (47)  into  Equation(45)  and  simplification  reveals 
that  the  covariance  matrix  is  the  inverse  of  the  Gauss-Newton  approximation  to  the 
Hessian, 


Vm  =  (JTVd  1J)~1  =  H-1  (48) 

The  above  equation  can  be  used  to  estimate  uncertainties  in  each  of  the  model  param¬ 
eters  and  can  also  be  used  to  define  confidence  ellipses  of  multiple  parameters.  The 
uncertainties  obtained  through  this  equation  are  linearized  estimates  of  the  true  un¬ 
certainties.  Before  considering  them  further  we  investigate  how  closely  these  linearized 
estimates  would  match  fully  non-linear  estimates.  If  the  minimum  of  the  objective 
function  is  m  then  linearized  estimates  are  equivalent  to  approximating  the  objective 
function  with 
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Depth  (m) 


Figure  16:  Objective  function  (solid  line)  and  Gauss-Newton  approximation  to  the  Hes¬ 
sian  (dashed  line)  about  the  inversion  solution  for  Target  575,  Fort  Ord. 


4>( m)  -  <t>{ rh)  «  ^(m  -  rh)TH(m  -  rh)  (49) 

If  the  above  equation  is  a  good  approximation  to  the  objective  function  about  the 
minimum  then  the  linearized  uncertainties  will  be  accurate  estimates  of  the  true  uncer¬ 
tainties.  We  have  compared  the  objective  functions  of  many  anomalies  and  usually  find 
a  situation  such  as  that  illustrated  in  Figure  (16).  The  linear  and  non-linear  estimates 
are  very  similar  for  parameters  related  to  the  dipole  moment,  and  less  so  for  parameters 
related  to  location  (especially  the  depth).  From  the  perspective  of  ordnance  identifica¬ 
tion  it  is  the  uncertainties  in  the  dipole  moment  that  axe  important,  so  we  can  conclude 
that  the  linearized  estimates  should  be  accurate  for  our  purposes.  We  note  in  passing 
that  including  the  second  order  terms  in  the  Hessian  has  very  little  impact  on  the  accu¬ 
racy  of  the  linear  approximation,  even  if  the  residuals  at  the  solution  are  large.  This  is 
readily  apparent  for  the  dipole  components  in  cartesian  coordinates  m  =  (rnx,my,mz) 
because  the  second  derivatives  with  respect  to  these  parameters  are  zero  (Equation  9). 

7.2  Confidence  ellipses  of  multiple  parameters 

We  are  principally  interested  in  the  joint-uncertainty  in  the  dipole  magnitude  and  an¬ 
gle  relative  to  the  Earth’s  field  because  these  two  parameters  control  our  decisions  on 
discrimination  and  identification.  A  good  method  for  visualizing  the  joint-uncertainties 
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is  to  calculate  confidence  ellipses;  these  are  contours  of  equal  confidence  in  the  model 
parameters.  They  are  usually  defined  to  encompass  a  certain  probability  that  the  model 
lies  within  the  defined  region.  For  instance,  a  one-sigma  (or  68.3%)  confidence  ellipse 
would  delineate  a  region  that  has  a  68.3%  chance  of  containing  the  true  model  parame¬ 
ters  (under  the  particular  statistical  assumptions  used  to  define  the  objective  function). 

To  define  the  confidence  ellipse  for  r  parameters  requires  that  the  uncertainties  in 
the  other  n  —  r  parameters  (assuming  there  are  n  parameters  in  total)  are  projected 
onto  the  r-dimensional  hyperplane  under  consideration.  In  turns  out  (Bard,  1974)  that 
this  projection  can  be  implicitly  achieved  by  calculating  the  covariance  matrix  through 
inversion  of  the  appropriate  r  x  r  subsection  of  the  Hessian  matrix.  For  example,  if  we 
are  interested  in  the  joint  confidence  ellipse  for  parameters  2  and  5  we  would  calculate 


Vw  = 

T  771 


#22  #25 

#52  #55 

The  principal  axes  of  the  ellipse  are  then  obtained 


J  (50) 

by  eigenvalue  decomposition  of  Vm, 


Vm  =  USUT  (51) 

where  S  is  a  diagonal  matrix  of  eigenvalues  and  U  is  an  orthonormal  matrix  of  eigen¬ 
vectors;  each  eigenvector  is  a  principal  axis  of  the  ellipse.  The  size  of  each  principal  axis 
is  given  by  the  reciprocal  of  the  corresponding  eigenvalue. 

The  final  step  is  to  find  the  ellipse  that  lies  on  the  boundary  of  a  specified  confidence 
region  (e.g.  68.3%  or  the  one-sigma  error).  If  the  sampling  distribution  is  normal,  unbi¬ 
ased  and  we  assume  that  the  covariance  matrix  Vm  is  known,  then  (m-m)TV“1(m— m) 
is  x2  distributed  with  r  degrees  of  freedom  (Bard,  1974).  The  boundary  of  the  ellipse 
corresponding  to  a  particular  confidence  level  is  found  by  calculating  the  value  of  the  x2 
distribution  with  r  degrees  of  freedom  that  corresponds  to  that  confidence.  For  example, 
with  2  degrees  of  freedom,  the  68.3%  confidence  limit  has  x2  ~  2.298.  The  contour  of 
(m  -  m)TV“1(m  -  m)  with  the  value  2.298  corresponds  to  the  boundary  of  the  68.3% 
confidence  region. 

Alternatively,  if  the  covariance  matrix  is  not  known  and  we  assume  that  all  observa¬ 
tions  are  independent  (so  that  Vd  =  a2 1)  then  the  confidence  level  is  obtained  through 
the  .Fr)n~r  distribution  (Bard,  1974).  In  that  case  the  value  of  x2  calculated  above  would 
be  replaced  with 


ra2Frtn-r{  7)  (52) 

where  7  is  the  desired  confidence  limit  (7  =  0.683  for  the  68.3%  confidence  level)  and  a2 
is  the  sample  variance.  If  the  optimum  model  returned  by  the  minimization  routine  is  m 
and  there  are  K  model  parameters  then  the  sample  variance  is  obtained  by  calculating 

a2  =  N^K^lFi^-d^]2  (53) 

i— 1 

Finally,  if  the  assumption  of  Gaussian  statistics  does  not  hold,  conservative  confidence 
limits  can  be  calculated  by  recourse  to  the  Bienayme-Chebyshev  inequality  (Bard,  1974). 
For  the  7  confidence  limit  this  replaces  the  value  of  x2  above  with 
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Figure  17:  Distribution  of  residuals  for  the  660  anomalies  at  Guthrie  Rd  with  correlation 
coefficient  greater  than  0.8. 


7.3  The  Guthrie  Road,  Montana  dataset 

The  distribution  of  residuals  for  the  660  anomalies  at  Guthrie  Road  with  correlation 
coefficients  greater  than  0.8  are  shown  in  Figure  (17).  The  curve  have  been  normalized 
so  that  they  have  an  integral  of  one;  thus  we  are  approximating  the  posterior  probability 
density  function  (posterior  PDF)  of  the  residuals.  The  best  fitting  Gaussian  distribution 
(with  a  =  0.99  nT)  provides  a  poor  fit  to  the  residuals.  In  particular,  large  residuals  are 
much  more  likely  than  predicted  by  a  normal  distribution.  In  Figure  (17)  we  also  show 
the  best  fitting  Ekblom  distribution  (Ekblom,  1973)  which  has  a  probability  density 
function  equal  to 


Pr(x)  =  cexp 


(x2  +  e2)p/2 
ap 


(55) 


where  c  is  a  constant  that  ensures  the  PDF  integrates  to  unity.  This  is  a  perturbed 
version  of  an  Lp-norm  and  we  found  best  fitting  values  of  a  =  1.16  x  10“4  nT,  e  = 
0.609  nT  and  p  =  0.224.  The  correspondence  between  the  observed  distribution  of 
residuals  and  the  Ekblom  PDF  is  excellent.  We  will  consider  the  good  fit  of  the  Ekblom 
PDF  further  in  the  discussion  section  of  this  report. 

As  the  residuals  are  clearly  non-Gaussian,  we  must  be  cautious  when  interpreting 
confidence  ellipses  calculated  under  the  normality  assumption.  The  shape  and  relative 
size  of  the  ellipses  are  likely  to  be  correct  but  their  absolute  values  are  probably  incorrect 
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Figure  18:  One  sigma  (68.3%)  and  two  sigma  (95.4%)  confidence  ellipses,  assuming 
Gaussian  statistics,  about  the  recovered  dipole  moments  for  (a)  76  mm  projectiles;  and 
(b)  81  mm  mortars. 

(i.e.  the  68.3%  confidence  ellipse  won’t  correspond  to  a  confidence  of  68.3%).  However, 
the  absolute  confidence  ellipses  should  provide  a  rough  estimate  of  the  uncertainties 
in  the  parameters.  They  should  also  provide  a  reliable  means  of  comparison  of  the 
uncertainties  in  recovered  parameters  for  different  anomalies. 

We  calculated  confidence  ellipses  for  the  dipole  magnitude  and  angle  relative  to  the 
Earth’s  field,  for  the  anomalies  due  to  76  mm  projectiles  and  81  mm  mortars  at  Guthrie 
Road  (Figure  18).  The  data- weighting  matrix  was  the  identity  matrix,  with  a  sample 
standard  deviation  calculated  from  the  residuals  by  Equation  (53).  The  figures  show 
the  one-sigma  (68.3%)  and  two-sigma  (95.4%)  confidence  ellipses  about  each  recovered 
dipole  moment  (assuming  Gaussian  statistics).  For  the  76  mm  projectiles  all  except  one 
anomaly  have  very  small  uncertainties.  The  one  outlier  was  incorrectly  identified  as  a 
105  mm  projectile.  For  the  81  mm  mortars  all  anomalies  have  very  small  uncertainties. 

We  also  calculated  confidence  ellipse  using  the  very  conservative  Bienayme-Chebyshev 
inequality  (Figure  19).  The  confidence  ellipses  are  significantly  larger  than  those  re¬ 
turned  under  the  normal  assumption.  If  the  real  confidence  ellipses  were  this  large,  then 
uncertainty  in  the  recovered  moment  would  be  a  contributing  factor  in  misidentifications 
(especially  for  76  mm  projectiles).  However,  the  real  confidence  ellipses  would  have  a 
size  intermediate  between  the  Gaussian  and  Bienayme-Chebyshev  estimates.  This  would 
imply  that  misidentification  is  sometimes,  though  probably  not  that  often,  caused  by 
uncertainty  in  the  recovered  moment.  Misidentifications  are  more  likely  to  occur  due  to 
modelling  deficiencies  and/or  remnant  magnetization. 

We  also  calculated  confidence  ellipses  (under  the  normal  assumption)  for  the  anoma¬ 
lies  due  to  metallic  debris  and  geology  (Figure  20).  For  the  metallic  debris  many  of  the 
uncertainties  are  quite  low  but  there  are  a  few  items  with  much  higher  uncertainties. 
For  the  anomalies  due  to  geology  the  uncertainties  are  generally  much  higher.  In  many 
cases  the  recovered  angle  is  around  60°  but  the  uncertainty  is  high.  This  is  a  critical  area 
for  making  discrimination  decisions  so  we  can  conclude  that,  for  geological  anomalies, 
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Figure  19:  Worst  case  scenario  for  one  sigma  (68.3%)  and  two  sigma  (95.4%)  confidence 
ellipses  about  the  recovered  dipole  moments  for  (a)  76  mm  projectiles;  and  (b)  81  mm 
mortars. 

uncertainty  in  dipole  angle  can  have  a  significant  impact  on  discrimination. 

8  DISCUSSION 

In  this  report  we  have  described  a  successful  inversion  method  for  recovering  the  po¬ 
sition  and  dipole  moment  from  a  total-field  magnetic  anomaly.  The  inversion  uses 
a  least-squares  formulation  and  an  interior-reflective  Newton  method  to  minimize  the 
objective  function.  A  starting  model  is  defined  using  an  analytical  method  based  on 
approximating  the  vertical  component  of  the  magnetic  field.  The  inversion  method  is 
very  fast  and  reliable;  our  Matlab  implementation  requires  only  a  few  seconds  to  recover 
the  parameters  for  each  anomaly. 

The  recovered  dipole  moments  can  be  used  to  make  discrimination  decisions  by 
recourse  to  the  phenomena  of  shock  demagnetization.  Complete  shock  demagnetization 
may  not  occur  for  all  ordnance  items  all  the  time.  Additionally,  non-ordnance,  even 
though  they  may  be  remnantly  magnetized,  might  be  orientated  so  that  their  dipole 
moments  lie  close  to  the  Earth’s  field.  This  means  that  it  is  not  possible  to  make  a 
definitive  ordnance/non-ordnance  decision.  Rather,  we  prefer  to  rank  anomalies,  with 
items  most  likely  to  be  UXO  higher  up  the  list.  We  investigated  three  such  ranking 
procedures;  (i)  dipole  angle,  (ii)  dipole  angle  and  magnitude  and  (iii)  percentage  of 
remnant  magnetism,  and  found  that  the  remnant  magnetization  ranking  was  the  best. 
With  this  method  all  ordnance  would  have  been  recovered  after  only  half  of  the  anomalies 
had  been  excavated.  This  included  85  anomalies  with  poor  dipole  fits  that  only  yielded 
one  UXO.  The  criterion  used  to  select  potential  UXO  may  need  to  be  revised  to  minimis 
the  number  of  anomalies  where  we  are  unable  to  make  reliable  inferences  about  the 
bodies  origin. 

Linearized  error  analysis  indicates  that  the  angle  between  the  recovered  dipole  mo¬ 
ment  and  the  Earth’s  field  for  ordnance  items  is  relatively  well  constrained.  Therefore, 
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Figure  20:  One  sigma  (68.3%)  and  two  sigma  (95.4%)  confidence  ellipses  about  the 
recovered  dipole  moments  for  (a)  metallic  debris;  and  (b)  geological  anomalies. 

uncertainty  in  the  recovered  dipole  moment  is  unlikely  to  strongly  influence  the  dis¬ 
crimination  ranking.  For  items  identified  as  geology,  the  confidence  ellipses  were  much 
larger.  Consequently,  the  position  of  these  items  in  the  ranking  is  strongly  influenced 
by  uncertainty. 

The  discrimination  method  was  very  successful  at  Guthrie  Road  but  could  not  be 
reliably  applied  to  the  seeded  sites  at  Fort  Ord.  This  is  a  consequence  of  remnant  mag¬ 
netization  and  emphasizes  the  importance  of  shock  demagnetization  to  the  viability  of 
the  method.  Additional  work  needs  to  be  undertaken  to  determine  the  circumstances 
where  shock  demagnetization  occurs.  For  instance,  the  recovered  dipole  moments  of 
several  of  the  81  mm  mortars  at  Guthrie  Road  had  relatively  large  dipole  angles,  in¬ 
dicating  the  presence  of  some  remnant  magnetism.  Perhaps  the  lower  impact  velocity 
of  mortars  compared  to  artillery  projectiles  causes  only  partial  shock  demagnetization; 
this  needs  to  be  investigated  further. 

Identification  of  the  ordnance  type  using  the  recovered  dipole  moment  is  difficult  due 
to  non-uniqueness.  The  dipole  moment  is  the  product  of  volume  with  the  magnetization 
vector.  When  the  volume  changes,  self-demagnetization  as  the  orientation  of  the  body 
is  varied  can  change  the  magnetization  vector  and  thus  produce  an  identical  or  similar 
dipole  moment.  However,  recognizing  that  there  are  only  a  finite  number  of  ordnance 
types  in  an  given  area,  makes  identification  using  the  recovered  dipole  moment  feasible. 
At  Guthrie  Road  correct  identification  could  be  achieved  about  55%  of  the  time. 

One  of  the  key  factors  that  influences  successful  identification  is  obtaining  the  correct 
trajectory  of  moment  magnitude  and  angle  for  a  given  ordnance  item.  We  use  a  spheroid 
approximation  to  calculate  this  trajectory.  Modelling  of  shapes  that  more  closely  ap¬ 
proximate  real  ordnance  should  be  undertaken  in  order  to  determine  the  reliability  of 
the  spheroid  approximation. 

A  significant  cause  of  misidentification  is  the  inability  of  the  dipole  moment  to  con¬ 
strain  the  orientation  of  the  causative  body.  In  an  attempt  to  overcome  this  problem  we 
developed  inversion  routines  that  recover  information  on  the  quadrupole  and  octupole 
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components  of  the  body.  The  dipole/octupole  method  we  developed  uses  a  spheroid 
approximation  of  ordnance  to  estimate  both  the  dipole  and  octupole  moments.  No  sig¬ 
nificant  improvement  in  identification  over  the  dipole  method  was  found.  This  could 
be  due  to  one  or  more  of  the  following:  (i)  remnant  magnetization;  (ii)  rapid  falloff  of 
the  octupole  component  with  distance  so  that  it  is  unable  to  constrain  the  orientation; 
and  (iii)  the  spheroid  provides  a  poor  approximation  to  the  octupole  moment  of  real 
ordnance. 

The  quadrupole  method  relaxed  the  requirement  of  top-bottom  symmetry  enforced 
by  a  spheroid  approximation  and  attempted  to  recover  the  dipole  and  quadrupole  mo¬ 
ments  of  an  arbitrary  axially  symmetric  body.  The  fits  to  the  data  were  improved 
slightly  but  no  advantages  for  identification  were  realized.  The  cause  is  likely  due  to  one 
or  more  of  the  following:  (i)  remnant  magnetization;  (ii)  the  inability  of  the  quadrupole 
component  to  constrain  the  orientation;  and  (iii)  inadequacy  of  the  model. 

The  failure  of  the  quadrupole  and  octupole  methods  to  improve  identification  may 
imply  that  analysis  of  magnetic  data  alone  does  not  provide  enough  information  to 
constrain  the  objects  orientation.  One  possible  method  to  provide  this  constraint  is  by 
joint  inversion  of  magnetic  and  electromagnetic  data;  either  single  or  multiple  channel 
time  domain  EM  or  frequency  domain  EM.  Most  EM  sensors  illuminate  the  target  from 
multiple  directions  and  are  consequently  more  sensitive  to  orientation  (in  contrast  the 
direction  of  the  Earth’s  field  is  fixed).  In  return,  the  magnetics  may  provide  valuable 
depth  constraints  that  are  difficult  to  obtain  with  EM  alone  (Pasion  and  Oldenburg 
2001). 

Linearized  error  analysis  was  used  to  estimate  individual  parameter  uncertainties  and 
joint  confidence  ellipses.  The  linear  approximation  provides  an  accurate  model  of  the 
objective  function  about  the  minimum.  However,  the  residuals  do  not  follow  a  Gaussian 
distribution  so  that  reliable  confidence  ellipses  could  not  be  defined.  Additionally,  the 
upper  bound  provided  by  the  Bienayme-Chebyshev  inequality  is  very  conservative  and 
probably  returns  confidence  ellipses  that  are  much  larger  than  required.  A  recommen¬ 
dation  for  future  work  is  to  attempt  to  calculate  confidence  ellipses  using  non-Gaussian 
statistics. 

The  residuals  were  found  to  follow  an  Ekblom  distribution,  yet  we  minimized  a 
measure  derived  from  Gaussian  statistics.  A  reformulation  of  the  problem  using  an 
Ekblom  distribution  as  a  measure  of  data  misfit  is  therefore  recommended.  The  resulting 
objective  function  would  need  to  be  solved  by  a  different  optimization  algorithm  than 
the  one  used  here.  Algorithms  for  the  solution  of  this  optimization  problem  are  well 
known  and  include  iteratively  reweighted  least  squares  and  Newton’s  method  (Watson, 
1980;  Ekblom,  1987).  Implementation  of  the  new  optimization  method  would  enable  a 
thorough  investigation  to  be  made  on  the  influence  the  choice  of  the  measure  of  data 
misfit  has  on  the  recovered  model.  This  is  an  issue  often  glossed  over  when  solving 
inverse  problems  and  this  report  was  no  exception. 

9  CONCLUSIONS 

A  dipole  inversion  method  was  developed  and  successfully  applied  to  the  discrimination 
and  identification  of  magnetic  anomalies  at  Guthrie  Road  Montana.  The  discrimination 
method  had  the  potential  to  reduce  the  number  of  excavations  by  half,  yet  still  yielded 
all  the  ordnance.  Identification  of  the  causes  of  the  anomalies  is  difficult  using  the 
recovered  dipole  moment  due  to  non-uniqueness.  However,  by  recognizing  that  there 
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are  a  finite  number  of  ordnance  types,  correct  identification  could  be  achieved  in  about 
55%  of  cases.  Extension  of  the  inversion  method  to  recover  either  quadrupole  of  octupole 
moments  did  not  improve  our  ability  to  discriminate  or  identify. 

We  recommend  the  following  areas  as  directions  for  future  research. 

1.  Shock  demagnetization:  The  extent  to  which  this  phenomena  occurs  should  be 
investigated  further  because  shock  demagnetization  is  central  to  the  success  of 
discrimination  using  magnetic  data. 

2.  Magnetic  modelling  of  ordnance:  The  trajectory  of  the  modelled  dipole  moment 
magnitude  and  angle  determines  our  identification  decisions.  Therefore,  the  agree¬ 
ment  between  the  moments  returned  by  a  spheroid  approximation  and  a  more 
realistic  ordnance  shape  should  be  investigated  further. 

3.  Joint  inversion  of  magnetics  and  electromagnetics:  Identification  difficulties  in 
magnetics  are  caused  by  an  inability  to  constrain  the  orientation  of  the  causative 
body.  The  ability  of  time  or  frequency  domain  electromagnetics  to  constrain  the 
orientation  needs  to  be  investigated.  Conversely,  magnetics  may  provide  valuable 
depth  information  that  can  be  exploited  by  the  EM  inversion. 

4.  Error  analysis  using  non-Gaussian  statistics:  The  residuals  do  not  follow  a  Gaus¬ 
sian  distribution  so  uncertainty  estimates  based  on  a  more  realistic  statistical 
distribution  should  be  investigated. 

5.  Robust  measures  of  data  misfit:  The  problem  should  be  reformulated  using  a  more 
realistic  measure  of  data  misfit  (the  Ekblom  norm),  and  methods  implemented  to 
solve  the  modified  inversion  problem. 
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A  MULTIPOLE  EXPANSION 


In  this  appendix  we  give  a  brief  derivation  of  the  multipole  expansion  of  a  magnetic 
body,  and  in  particular  derive  the  equations  relevant  to  a  spheroid.  Our  derivations 
closely  follow  those  given  in  McFee  (1989). 

The  scalar  potential,  <fi(r)  of  a  body  with  magnetization  M  is  given  by 

#w"sXM'v(?)‘n'  (A-1) 

where  r*  is  the  distance  from  the  source  to  the  observation  point  (Figure  1).  Expanding 
the  gradient  and  using  the  divergence  theorem  we  find 

m  =  h{Js?MdS-Jv?v' MdV)  (A-2) 

If  the  induced  field  inside  the  body  is  uniform  and  parallel  (which  it  is  for  a  spheroid), 
which  eliminates  the  second  term  in  the  above  equation  because  then  V  •  M  =  0.  The 
multipole  method  proceeds  by  expanding  p  as  a  Taylor  series  about  the  origin, 
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where  repeated  indices  imply  that  summation  convention  is  to  be  used.  The  scalar 
potential  from  a  2n  pole  with  moment  is  given  by  the  expression 
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from  which  it  follows  that  the  first  four  moments  are 
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In  order  the  moments  are  the  scalar  monopole  m^°\  the  vector  dipole  m^l\  the  rank  2 
tensor  quadrupole  and  the  rank  3  tensor  octupole  m^3\  Magnetic  monopoles  do 
not  exist  which  means  that  the  monopole  term  is  zero.  Expansion  of  the  dipole  term 
reveals  that  the  dipole  moment  is  the  product  of  the  magnetization  with  the  volume, 
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m(i)  _  MV  (A-9) 

Due  to  the  symmetry  of  the  spheroid  the  quadrupole  term  is  zero.  After  the  dipole,  the 
next  non-zero  term  is  the  octupole;  specification  of  the  octupole  components  is  tedious 
so  we  refer  the  interested  reader  to  McFee  (1989). 

The  magnetic  field  (in  the  i-th  direction)  due  to  the  dipole  component  is  given  by 
the  expression 


(A- 10) 


Although  zero  for  the  spheroid  we  will  use  the  quadrupole  term  later,  its  field  is, 


Bi2)  =  ^5  +  mfi)  +  5r  2XiXjXkmf^  (A-ll) 
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and  for  the  octupole  the  field  is, 


Bi3)  =  g^5  (3mS  -  15r  2  [xiximfkk  +  xJxkm<§l}  +  35r  4xixjxkxlmf^  (A- 12) 

It  remains  to  find  the  magnetization  vector  for  a  spheriod.  A  solution  of  the  boundary 
value  problem  (Stratton,  1941)  shows  that  the  demagnetization  factors  are 


F  =  Mr-1 
1  1  4-  Oii(fir  —  l)/2 

where  /zr  is  the  relative  permeability  (/z  =  /zr/z0), 


e(e  +  E) 
ai=  a2  =  — 5 — — 
el  -  1 

and 


(A-13) 


(A-14) 


with 


— 2e(e  ^  ~h  E} 
e2  -  1 


(A-15) 


log(e  -  %/e2  -  1) 

E  =  VS-i 

for  a  prolate  spheroid  (e  >  1)  and 

arctan(-7=j)  —  tt/2 

fcj  =  - ,  - 

%/T^ 

for  an  oblate  spheroid  (e  <  1).  For  the  special  case  of  a  sphere  (e  =  1)  then 


(A- 16) 


(A- 17) 


2 

=  a2  =  a3  =  - 


(A- 18) 


The  demagnetization  factors  together  with  the  Earth’s  field  B0  (in  spheroid  centered 
coordinates),  determine  the  strength  of  the  induced  magnetization, 


M  =  ^ FB0 


(A- 19) 


where  we  have  written  F  as  a  3  x  3  diagonal  matrix  with  {F2)F2,Fy)  along  the  diagonal 
(because  F\  =  F2). 
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A.l  Multipole  expansion  for  an  axi-symmetric  body 

In  going  from  the  general  equation  (A-2)  for  a  magnetized  body  to  that  for  a  spheroid, 
we  used  the  fact  that  the  magnetization  inside  the  body  was  constant  so  that  V  M  =  0. 
For  a  general  axi-symmetric  body  this  is  no  longer  true  and  there  is  an  extra  contribution 
to  each  of  the  moments.  For  the  dipole, 

mg)  =  J  •  dS  -  J  Xi  V  •  MdV  (A-20) 

and  using  the  identity 

J  4>M  ■  V<fidV  =  J  </>M  •  dS  -  J  <t>y  -MdV,  (A-21) 

we  find 

mp)  =  [  M-VxidV 

Jv 

=  Jv  MidV  (A-22) 

Thus  the  dipole  moment  is  just  the  integration  of  the  magnetization  over  the  volume 
of  the  body.  Most  ordnance  items  are  either  cylinder  like  with  a  pointed  front  (e.g. 
most  artillery  projectiles)  or  close  to  a  spheroid  shape  (e.g.  mortars,  the  fins  can  be 
ignored  as  they  are  usually  aluminum).  For  an  infinite  cylinder  the  field  lines  are  con¬ 
stant  within  the  body.  Therefore,  for  a  projectile  any  changes  in  magnetization  will  be 
concentrated  at  the  ends  of  the  body,  and  will  generally  make  only  a  small  contribution 
to  the  multipole  moments.  For  a  mortar,  the  shape  is  close  to  second  order  (Stratton, 
1941)  and  hence  the  magnetization  direction  will  be  approximately  constant.  Thus, 
changes  in  magnetization  direction  should  not  make  a  significant  contribution  to  the 
observed  magnetic  anomaly. 

Using  an  analogous  development  to  the  dipole,  the  quadrupole  moment  can  be  shown 
to  be  equal  to 

mg)  =  M  •  V(xlXj)dV 

=  J  ( MiXj  +  MjXi)dV  (A-23) 

Assuming  that  changes  in  magnetization  can  be  ignored  implies  that 

=  Mi  [  dxs  f  XjdS  +  Mj  [  dx 3  f  XidS  (A-24) 

3  Jl  Js(x  3)  Jl  Js{xz) 

where  L  is  the  length  of  the  body  parallel  to  the  symmetry  axis,  z  —  x 3  and  S(x 3)  is  a 
slice  of  the  body  perpendicular  to  the  axis  of  symmetry.  Due  to  symmetry  each  of  the 
integrals  over  £(23)  are  zero  for  i  or  j  =  1  or  2,  from  which  it  follows  that 

mff  =  0  for  any  i,j  =  1  or  2  (A-25) 

(2)  (2) 

The  quadrupole  moment  is  symmetric  m\j  —  rrv^  so  that  the  full  moment  tensor  can 
then  be  written  as  (where  we  have  dropped  the  (2)  subscript), 


B  CONSTRAINTS  ON  DIPOLE  ORIENTATION 

The  angle  the  induced  dipole  moment  in  a  spheroid  makes  with  the  Earth’s  field  cannot 
exceed  a  certain  value.  In  this  appendix  we  derive  the  equations  which  specify  this  upper 
bound.  We  assume  the  Earth’s  field  is  at  an  angle  8  to  the  spheroid  symmetry  axis  and 
align  the  z-axis  of  our  coordinate  system  with  the  symmetry  axis.  We  can  then  choose 
the  x-axis  so  that  it  is  perpendicular  to  the  Earth’s  field.  In  this  coordinate  system 

bD  =  6o(0,  sin  0,  cos  0)T  (B-l) 

The  spheroid  diameter  a  and  eccentricity  e  determine  the  demagnetization  factors,  F2 
and  F3,  so  that  the  induced  dipole  moment  m  is, 

m  =  V  b0( 0,  F2  sin  8 ,  F$  cos  0)T  (B-2) 

where  V  is  the  volume  of  the  spheroid.  The  angle  ip  between  the  Earth’s  field  and  the 
dipole  moment  is 


cos  ip  = 


m 


llbolllHI 

which,  on  substitution  of  Equations  (B-l)  and  (B-2),  becomes 

F2  sin2  8  +  F%  cos2  8 


cos  ip  = 


(B-3) 


(B-4) 


After  differentiating  this  equation  by  8 ,  setting  the  result  to  zero  and  extensive 
algebraic  manipulation  we  find 


8  —  arctan 


(  I Fj  +  FsFj-2Fhjf\ 
F$  +  F2Fi-2F3Fi) 


(B-5) 


This  orientation  of  the  spheroid  symmetry  axis  relative  to  the  Earth’s  field  makes  the 
angle  between  the  induced  dipole  and  the  Earth’s  field  a  maximum.  This  maximum 
angle  is  obtained  by  substituting  Equation  (B-5)  into  Equation  (B-4). 


C  AMBIGUITY  IN  THE  DIPOLE  SOLUTION 

Once  the  sensor  distance  exceeds  a  few  body  lengths  distance  away  from  the  spheroid, 
the  field  is  essentially  dipolar.  In  such  cases,  inversion  for  spheroid  dimensions  has  an 
inherent  ambiguity:  many  different  spheroid  diameters  and  aspect-ratios  can  reproduce 
the  observed  dipole.  To  characterize  this  ambiguity  we  chose  our  coordinate  system  so 
that  the  z-axis  is  aligned  with  the  Earth’s  field, 

B0  =  (0, 0,  B0)t  (C-l) 

We  want  to  determine  if  the  dipole  m  from  a  particular  spheroid  (diameter  a  and 
aspect-ratio  e)  can  match  an  observed  dipole  moment  m.  Now  the  induced  dipole  for 
the  spheroid  is  given  by  the  following  equation 

m  =  VArFAB0  (C-2) 
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v  * 


Expanding  this  equation  we  find 


m  = 


Bo'nea 3 
6fi0 


(F2  —  Fs)  cos  9  sin  6  sin  xp 
(F2  -  Fs)  cos  9  sin  9  cos  xp 
F2  cos2  9  4-  F$  sin2  9 


(C-3) 


Now  we  want  to  find  a,  e,  9  and  xp  so  that  m  =  m.  We  can  find  an  expression  for  the 
spheroid  diameter  by  equating  the  norms  of  the  two  dipole  moments,  resulting  in 


For  xp  we  use, 


which  leads  to 


_ fyolMI _ 

B0ne  ^FicosH  +  Fi  sin20 


rxi\  mi  sin  xp 

—  =  t—  = - -  =  tan  xp 

m2  m2  cos  xp 


(C-4) 


(0-5) 


xp  =  arctan 


(C-6) 


This  implies  that  the  spheroid  has  the  same  azimuth  as  the  dipole  moment. 

Determination  of  9  is  a  little  more  difficult.  We  start  by  using  the  second  row  of 
Equation  (C-3)  and  substitute  in  Equation  (C-4)  for  the  diameter, 


m2  _  (F2  -F3)cosflsin6> 

cos^IMI  ^F22cos2£  +  F32sin20 

Note  that  if  cos  xp  =  0  we  would  use  the  first  row  of  Equation  (C-3)  to  derive  an  equation 
in  terms  of  mi  and  sin^  ■  Let  the  LHS  of  Equation  (C-7)  be  cj,  expand  out  the  square- 
root  on  the  denominator  of  the  RHS,  and  use  the  substitution 


sin2  9  —  1  -  cos2  9  (C-8) 

After  some  algebraic  manipulation  Equation  (C-7)  reduces  to 

(F2  -  F3)2  cos4  e  p(Ff  -  F32)  -  (F2  -  F3)2]  cos2  e  +  u>2F$  (C-9) 

This  is  a  quadratic  equation  in  cos2  9  which  has  a  real  solution  if 

A  ~b2-  4ac  >  0  (C-10) 

where 


a  =  (f2  -  f3)2  (c-ii) 

b  =  w\F2  -  F32)  -  (F2  -  F3)2  (C-12) 

c  =  a ;2F32  (C-13) 


The  solution  in  terms  of  6  is  then 
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+*  * 


0  =  zb  arccos 


-b±VE 

2a 


(C-14) 


The  =L  terms  in  the  above  equation  emphasize  that  when  A  >  0  there  are  actually  two 
possible  solutions. 

Determination  of  which  spheroids  can  reproduce  the  observed  dipole  proceeds  as 
follows.  Firstly,  we  chose  an  aspect-ratio  e  and  calculate  the  demagnetization  factors  in 
the  usual  way.  We  then  substitute  Equations  (C-ll)  to  (C-13)  into  Equation  (C-10)  to 
determine  if  a  solution  can  be  found  for  this  aspect-ratio.  If  it  can,  we  calculate  6  by 
Equation  (C-14),  and  finally  we  calculate  the  spheroid  diameters  (usually  there  will  be 
two)  by  Equation  (C-4). 


Diameter  (mm) 

Aspect  ratio 

Angle  (degrees) 

A 

82 

5 

40.5 

B 

138 

5 

84.0 

C 

138 

2.82 

57.2 

D 

327 

0.06 

53.4 

E 

327 

0.31 

33.4 

Table  1:  Spheroid  dimensions,  and  their  angles  relative  to  the  Earth’s  field,  that  produce 
the  same  dipole  moment  as  a  105  mm  projectile  at  45°  inclination. 


76  mm 

81  mm 

Large  pieces 

Schrapnel 

Metal 

Geology 

Total 

All  items 

34 

49 

19 

217 

377 

126 

822 

Failed  fits 

0 

1 

1 

13 

41 

29 

85 

Table  2:  Number  of  items  found  during  excavation,  and  the  number  of  fits  with  c  <  0.7 


Angle 

Dig 

Leave 

Ordnance  found 

Ordnance  missed 

False  alarms 

55 

449 

373 

78 

5 

371 

60 

485 

337 

81 

2 

404 

65 

522 

300 

82 

1 

440 

70 

553 

269 

82 

1 

471 

75 

605 

217 

83 

0 

522 

Table  3:  Guthrie  Road  discrimination  results  for  intact  ordnance  items.  Discrimination 
used  the  dipole  angle  only. 
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Angle 

Dig 

Leave 

Ordnance  found 

Ordnance  missed 

False  alarms 

55 

449 

373 

91 

11 

358 

60 

485 

337 

99 

3 

386 

65 

522 

300 

101 

1 

421 

70 

553 

269 

101 

1 

452 

75 

605 

217 

102 

0 

503 

Table  4:  Guthrie  Road  discrimination  results,  including  large  pieces  of  ordnance  items. 
Discrimination  used  the  dipole  angle  only. 


Mag.  (A/m) 

Dig 

Leave 

Ordnance  found 

Ordnance  missed 

False  alarms 

0.06 

383 

439 

81 

2 

302 

0.05 

443 

379 

83 

0 

360 

Table  5:  Guthrie  Road  discrimination  results  for  intact  ordnance  items.  Discrimination 
used  the  dipole  angle  (75°)  and  magnitude. 


Mag.  (Am2) 

Dig 

Leave 

Ordnance  found 

Ordnance  missed 

False  alarms 

0.06 

383 

439 

90 

12 

293 

0.05 

443 

379 

95 

7 

348 

0.04 

493 

329 

98 

4 

395 

0.03 

558 

264 

100 

2 

458 

0.02 

597 

1  225 

102 

0 

495 

Table  6:  Guthrie  Road  discrimination  results,  including  large  pieces  of  ordnance  items. 
Discrimination  used  the  dipole  angle  (75°)  and  magnitude. 


Remnant  (%) 

Dig 

Leave 

Ordnance  found 

Ordnance  missed 

False  alarms 

10 

228 

594 

46 

37 

182 

20 

306 

516 

71 

12 

235 

30 

350 

472 

79 

4 

271 

40 

384 

438 

82 

1 

302 

50 

415 

407 

83 

0 

332 

Table  7:  Guthrie  Road  discrimination  results  for  intact  ordnance  items.  Discrimination 
used  the  dipole  angle  (75°),  magnitude  (<  0.05  Am2)  and  remnant  magnetism. 


Item 

Diameter  (mm) 

Eccentricity 

60  mm  mortar 

60 

3.78 

76  mm  projectile 

76 

3.57 

81  mm  mortar 

81 

3.54 

90  mm  projectile 

90 

5.0 

105  mm  projectile 

105 

3.71 

155  mm  projectile 

155 

4.0 

Table  8:  Ordnance  sizes  used  for  the  Guthrie  Road  identification  results. 
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% 


Number 

Item 

%a  (m) 

Zr  (m) 

c 

Angle 

Interp 

Rank 

620 

Empty  Hole 

0.00 

-0.47 

0.79 

127.4 

Scrap 

617 

Empty  Hole 

0.00 

-0.29 

0.75 

97.0 

Scrap 

599 

Fragment 

0.30 

-0.92 

0.88 

86.1 

Scrap 

598 

Fragment 

0.30 

-0.62 

0.90 

114.4 

Scrap 

597 

Fragment 

0.15 

-0.28 

0.97 

48.8 

Signals 

596 

Fragment 

0.30 

-0.41 

0.91 

55.2 

22  mm 

616 

35  mm 

0.61 

-0.56 

0.93 

74.0 

Scrap 

608 

35  mm 

0.46 

-0.23 

0.85 

59.0 

22  mm 

2 

607 

35  mm 

0.30 

-0.23 

0.82 

150.1 

Scrap 

606 

35mm 

0.15 

-0.25 

0.98 

33.6 

37  mm 

6 

605 

35  mm 

0.46 

-0.95 

0.91 

52.4 

Signals 

3 

604 

35  mm 

0.61 

-1.50 

0.81 

65.3 

Scrap 

603 

35  mm 

0.15 

-0.37 

0.90 

55.2 

22  mm 

2 

600 

35  mm 

0.30 

-0.43 

0.86 

68.8 

Scrap 

593 

35  mm 

0.61 

-0.53 

0.92 

38.0 

22  mm 

3 

592 

35  mm 

0.30 

-0.44 

0.97 

87.1 

Scrap 

591 

35  mm 

0.61 

-0.34 

0.95 

105.6 

Scrap 

590 

35  mm 

0.46 

-0.38 

0.96 

37.5 

60  mm  mortar 

2 

613 

2.36”  rocket 

0.61 

-0.65 

0.96 

57.1 

Rocket  2  .36” 

1 

602 

2.36”  rocket 

0.46 

-0.59 

0.98 

55.1 

Signals 

3 

594 

2.36”  rocket 

0.61 

-0.89 

0.96 

50.8 

Rocket  3  .5” 

4 

586 

2.36”  rocket 

0.61 

-2.79 

0.96 

118.2 

Scrap 

615 

Signals 

0.15 

-0.26 

0.97 

39.7 

37  mm 

2 

614 

Signals 

0.30 

-0.61 

0.82 

82.6 

Scrap 

611 

Signals 

0.15 

-0.33 

0.89 

43.3 

22  mm 

3 

589 

Signals 

0.30 

-1.08 

0.89 

20.1 

22  mm 

8 

588 

Signals 

0.30 

-1.51 

0.95 

31.6 

81  mm  mortar 

6 

610 

Hand  grenade 

0.15 

-0.08 

0.93 

49.9 

22  mm 

12 

609 

Hand  grenade 

0.30 

-0.58 

0.77 

38.6 

22  mm 

12 

587 

Hand  grenade 

0.30 

-0.88 

0.99 

121.5 

Scrap 

601 

Rifle  grenade 

0.46 

-1.33 

0.78 

77.9 

Scrap 

Table  9:  Fort  Ord  identification  results  for  Kl,  where  za  is  the  actual  burial  depth,  zr 
is  the  recovered  depth,  c  is  the  correlation  coefficient  between  observed  and  predicted 
data  and  angle  is  the  angle  between  the  Earth’s  field  and  the  recovered  dipole  moment. 
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**  %» 


Number 

Item 

Za  (m) 

Zr  (m) 

c 

Angle 

Interp 

Rank 

565 

Fragment 

-0.15 

-0.14 

0.93 

92.1 

Scrap 

566 

Fragment 

-0.30 

-0.34 

0.98 

113.5 

Scrap 

567 

Fragment 

-0.15 

-0.14 

0.95 

52.3 

22  mm 

568 

Fragment 

-0.30 

-0.28 

0.89 

37.3 

37  mm 

578 

Fragment 

-0.15 

-0.22 

0.92 

32.5 

60  mm  mortar 

579 

Fragment 

-0.30 

-1.82 

0.86 

90.1 

Scrap 

581 

Fragment 

-0.30 

-0.57 

0.90 

36.4 

22  mm 

580 

25  mm 

-0.61 

-0.59 

0.82 

72.1 

Scrap 

553 

35  mm 

-0.61 

-0.67 

0.95 

29.1 

35  mm 

1 

559 

35  mm 

-0.61 

-0.77 

0.89 

33.2 

75  mm 

8 

575 

35  mm 

-0.30 

-0.26 

0.97 

154.0 

Scrap 

563 

37  mm 

-0.46 

-0.20 

0.81 

55.6 

22  mm 

3 

583 

37  mm 

-0.61 

-0.93 

0.85 

47.7 

35  mm 

5 

584 

37  mm 

-0.61 

-1.03 

0.89 

61.8 

Scrap 

554 

2.36”  rocket 

-0.61 

-0.89 

0.97 

21.6 

90  mm 

9 

558 

2.36”  rocket 

-0.61 

-0.81 

0.94 

49.5 

105  mm 

3 

564 

75  mm 

-0.76 

!  -0.67 

0.96 

!  12.6 

75  mm 

1 

572 

75  mm 

-0.91 

-1.85 

0.86 

85.6 

Scrap 

560 

Signals 

-0.30 

-0.50 

0.81 

124.4 

Scrap 

576 

Signals 

-0.15 

-0.16 

0.95 

99.8 

Scrap 

573 

3”  Stokes 

-1.02 

-0.63 

0.86 

60.0 

35  mm 

10 

577 

3”  Stokes 

-1.02 

-1.22 

0.89 

19.2 

155  mm 

2 

585 

3”  Stokes 

-0.91 

-0.66 

0.87 

16.2 

35  mm 

9 

561 

81  mm 

-0.91 

-1.13 

0.96 

39.5 

81  mm  mortar 

1 

571 

81  mm 

-0.91 

-0.85 

0.94 

34.8 

90  mm 

2 

557 

90  mm 

-0.76 

-0.82 

0.97 

21.3 

90  mm 

1 

570 

90  mm 

-0.76 

-0.73 

0.97 

20.2 

90  mm 

1 

562 

105  mm 

-1.22 

-1.16 

0.94 

59.7 

Rocket  2.36” 

3 

582 

105  mm 

-1.22 

-1.50 

0.96 

17.8 

105  mm 

1 

552 

Hand  grenade 

-0.30 

-0.44 

0.92 

44.4 

Signals 

12 

Table  10:  Fort  Ord  identification  results  for  K2,  where  za  is  the  actual  burial  depth,  zr 
is  the  recovered  depth,  c  is  the  correlation  coefficient  between  observed  and  predicted 
data  and  angle  is  the  angle  between  the  Earth’s  field  and  the  recovered  dipole  moment 
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